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Chapter  1 

Executive  Summary 


Hyperspectral  Imagery  (HSI)  remains  the  most  advanced  tool  for  passive  remote  sensing. 
A  hyperspectral  data  cube  can  contain  a  vast  range  of  spatial,  spectral,  and  temporal 
information.  Interest  in  Hyperspectral  imagery  is  creating  markets  for  such  applications 
as  military  intelligence,  tactical  military  operations,  geology,  forestry,  cultural  studies,  and 
environmental  studies.  These  disparate  remote  sensing  communities  require  sophisticated 
hyperspectral  techniques  to  solve  intrinsically  difficult  problems  such  as  automated  local  area 
change  detection,  small  target  detection  in  clutter,  material  identification  at  the  sub-pixel 
scale,  target  class  discrimination  for  targeting,  and  geophysical  or  biochemical  parameter 
estimation  for  natural  resource  development  or  environmental  monitoring. 

However  there  are  limitations.  Hyperspectral  imaging  systems  are  usually  of  lower  spatial 
resolution  than  multi-spectral  or  panchromatic  (broadband)  imaging  systems.  The  low 
to  medium  spatial  resolution  of  hyperspectral  imagery  limits  its  application  to  real-world 
problems  such  as  scene  spatial  change  detection  or  image  analysis. 

In  completing  this  Phase  2  SBIR,  ATK  Mission  Research  with  help  from  the  University 
of  Dayton,  has  accomplished  the  intended  goal  of  the  SBIR  program  by  developing  a  com¬ 
mercially  marketable  Hyperspectral  Image  Enhancement  system.  This  system  implements 
a  coherent  and  rigorous  physics  model  for  HSI  spatial  resolution  enhancement  on  a  Mercury 
Multi-Computer  platform.  While  initial  development  was  accomplished  using  the  8-node 
Mercury  AdapDev  1280  computer  system,  we  have  demonstrated  that  our  code  easily  scales 
to  larger  vector  processing  clusters  such  as  the  96  (Titan)  and  64  (Hawk)  node  Mercury 
systems  residing  at  the  Wright  Patterson  Air  Force  Base  Computer  Facility. 

Our  SBIR  research  effort  has  led  to  the  completion  of  a  Maximum  A  Posteriori  (MAP) 
algorithm  for  constructing  a  high-spatial  resolution  hyperspectral  enhancement  from  a  low- 
spatial  resolution  hyperspectral  image,  by  incorporating  registered  higher-resolution  imagery 
from  an  auxiliary  sensor.  Although  we  have  focused  on  the  enhancement  of  a  hyperspectral 
image  using  high-resolution  panchromatic  data,  the  estimation  software  delivered  allows  for 
any  number  of  spectral  bands  to  be  employed  from  the  auxiliary  sensor. 
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The  technique  is  suitable  when  some  correlation  exists  between  the  auxiliary  image 
and  the  image  being  enhanced.  For  example,  LWIR  HSI  may  be  improved  using  MWIR 
broadband  imagery  as  shown  in  Figures  1.1  and  1.2.  Using  principle  component  analysis, 
we  have  shown  that  our  MAP  estimator  produces  a  hyperspectral  image  with  enhanced 
sub-pixel  content  evident  not  only  in  the  first  principal  component  image,  but  in  lower 
components  as  well.  Typical  principal  component  substitution  methods  seek  to  enhance  only 
the  first  principal  component  and  do  not  enhance  lower  components  at  all.  Any  enhancement 
of  these  low-signal  lower  components  indicates  a  promising  result. 


Figure  1.2:  SEBASS  imagery:  MWIR  broadband 


Figure  1.3:  SEBASS  imagery:  MAP  estimate  of  band  50  LWIR  using  MWIR  broadband 
imagery 
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Figure  1.4:  SEBASS  imagery:  True  band  50  LWIR 
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We  believe  that  our  MAP  approach  represents  the  best  results  from  the  most  theoreti¬ 
cally  sound  method  proposed  to  date  for  merging  hyperspectral  imagery  with  an  auxiliary 
sensor.  In  any  applications  where  panchromatic  sharpening  is  used  for  human  interpreta¬ 
tion,  our  method  offers  improved  performance.  A  technical  strength  of  our  MAP  estimation 
framework  is  that  it  allows  for  the  incorporation  of  alternate  methods  of  parameter  estima¬ 
tion  into  the  existing  code.  This  modifiable  estimation  framework  will  allow  AFRL/SN  to 
explore  new  methods  for  estimating  statistical  parameters  which  hold  the  promise  of  still 
better  image  enhancement. 

The  original  proposal  objectives  follow.  We  believe  that  we  have  successfully  addressed 
the  tasks  defined  in  the  initial  proposal  and  outlined  below.  As  the  tasks  evolved,  we 
reassigned  resources  to  Objective  1  in  order  to  meet  the  primary  emphasis  of  the  contract, 
to  develop  the  prototype  system. 

Objective  1:  Develop  a  prototype  near  real  time  hyperspectral  resolution  enhancement 
system  that  implements  the  algorithm  developed  during  Phase  I. 

Complete:  This  task  represented  80%  of  our  technical  effort.  The  Mercury  1280  AdapDev 
system  we  are  delivering  is  a  prototype  multi-computer  resolution  enhancement  system 
incorporating  the  final  and  improved  version  of  our  algorithm. 

Objective  2:  Further  develop  the  MAP  estimator  to  take  advantage  of  any  a  priori  spectral 
information  such  as  known  spectral  signatures  in  the  scene. 

Complete:  The  final  algorithm  incorporates  capabilities  for  incorporating  user-selectable 
multi-spectral  imagery  and  for  outputting  Principal  Component  imagery.  We  have  also 
included  the  capability  for  users  to  incorporate  their  own  linear  observation  models 
for  combinations  of  endmember  spectra. 

Objective  3:  Incorporate  techniques  for  Temperature/Emissivity  Separation  (TES)  for 
calibrated  high-resolution  hyperspectral  data.  The  goal  is  to  recover  land  surface 
temperature  from  the  calibrated  resolution  enhanced  hyperspectral  data. 

Complete:  Although  not  implemented  in  the  final  product,  a  technique  was  researched 
and  detailed  in  Appendix  C.  It  is  anticipated  that  future  customers  will  have  their 
own  favored  TES  algorithms  that  we  can  implement  as  required. 

Objective  4:  Incorporate  multi-frame  enhancement  techniques  to  be  applied  to  the  broad¬ 
band  imagery  prior  to  hyperspectral  fusion. 

Complete:  Again,  our  final  algorithm  includes  the  capability  for  users  to  incorporate  their 
own  linear  observation  models  for  combinations  of  endmember  spectra.  We  have  not 
specifically  included  multi-frame  enhancement  techniques  in  the  final  product  as  most 
hyperspectral  and  high-resolution  panchromatic  imaging  systems  incorporate  scanning 
technology  rather  than  focal  plane  arrays  suitable  for  multi- frame  enhancement. 
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Objective  5:  Incorporate  Air  Force  Research  Laboratories  atmospheric  radiance  and  trans¬ 
mission  models  to  allow  atmospheric  correction  for  raw  hyperspectral  data. 

Complete:  During  this  task  we  ran  a  NASA/JPL  parallelized  version  of  MODTRAN  3 
on  our  multi-node  Beowulf  Cluster.  We  stopped  short  of  fully  incorporating  parallel 
FLAASH  code  into  the  AdapDev  due  to  funding  constraints.  Our  implementation 
plan  was  to  first  operate  the  parallel  FLAASH  code  on  a  multi-node  Beowulf  cluster 
and  then  to  migrate  the  code  to  the  Mercury  platform.  We  still  believe  that  this  is  a 
feasible  approach  and  propose  to  investigate  it  as  a  follow-on  step. 

Objective  6:  Incorporate  image  registration  capability  in  the  prototype  to  align  low- resolution 
hyperspectral  imagery  with  high-resolution  broadband  imagery. 

Complete:  We  started  this  task  with  a  simple  registration  error  analysis.  Our  results 
demonstrated  that  if  the  panchromatic  auxiliary  image  is  perturbed  by  even  1  full 
pixel  from  the  HSI  data  our  MAP  algorithm  is  no  longer  any  more  effective  than  simple 
interpolation.  Even  a  registration  error  of  a  fraction  of  a  pixel  degrades  performance 
substantially.  Because  the  MAP  algorithm  is  so  sensitive  to  registration  we  conclude 
that  our  Image  Enhancement  product  is  best  suited  to  improving  HSI  collocated  with 
the  auxiliary  sensor  and  sharing  common  viewing  optics.  This  platform  architecture 
will  reduce  registration  error  upfront. 

Future  Work 

1.  Explore  other  forms  for  the  critical  MAP  algorithm  conditional  covariance  matrix 
and  imply  other  assumptions  regarding  the  nature  of  the  high-resolution  hyper-pixels. 
There  is  on-going  work  in  the  area  of  Blind  Source  Separation  that  may  help  with  the 
hyperspectral  pixel  unmixing  problem. 

2.  Develop  a  covariance  matrix  partitioning  scheme  that  increases  MAP  algorithm  per¬ 
formance,  subject  to  speed  and  memory  tradeoffs.  The  results  will  depend  on  the 
specific  scene  used  and  the  relative  quantities  of  the  spectra  present. 

3.  Incorporate  the  parallelized  version  of  the  atmospheric  correction  code  FLAASH  into 
the  MAP  enhancement  algorithm  on  a  Mercury  platform. 
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Chapter  2 


Introduction  and  Overview 


The  inherent  trade-off  between  spatial  and  spectral  resolution  in  hyperspectral  sensors  has 
prompted  the  development  of  remote  sensing  systems  that  include  low-resolution  hyperspec¬ 
tral  coupled  with  high- resolution  panchromatic  and/or  multispectral  imaging  subsystems. 
An  example  is  the  NASA  Earth  Observer  1  satellite,  which  includes  a  30  m  hyperspectral 
sensor  and  a  10  m  panchromatic  imager.  Commercial  panchromatic  satellite  imagery  ap¬ 
proaching  1  m  spatial  resolution  is  also  available.  This  provides  the  opportunity  to  jointly 
process  the  hyperspectral  and  higher  resolution  panchromatic  imagery  to  potentially  achieve 
improved  detection  and/or  classification  performance,  improved  change  detection,  and  im¬ 
proved  visual  imagery  for  human  interpretation. 

A  variety  of  techniques  have  been  presented  in  the  literature  for  merging  imagery  of  dif¬ 
ferent  spatial  and  spectral  resolution  [1-15].  Many  of  these  techniques  have  been  designed 
to  sharpen  multispectral  imagery  for  human  interpretation  using  broadband  panchromatic 
data.  Component  substitution  methods  transform  the  multispectral  imagery  and  replace 
one  component  with  the  broadband  high-resolution  imagery  [4,6,16].  Commonly  used  trans¬ 
formations  are  intensity- hue-saturation  (IHS),  where  the  intensity  is  replaced,  or  principal 
component  analysis  (PCA),  where  the  first  principal  component  is  replaced.  Clearly,  infor¬ 
mation  in  the  lower  components,  that  may  be  critical  in  classification  and  detection,  is  not 
enhanced  with  such  an  approach. 

High  pass  techniques  add  high  spatial  frequency  content  from  the  high-resolution  image 
to  the  bands  of  the  low-resolution  data  [3,5,8, 9, 16].  The  technique  in  [1]  uses  a  statistical 
approach  that  adds  a  linear  combination  of  the  high-resolution  data  to  a  pixel  replicated 
version  of  the  low-resolution  imagery.  This  technique  was  designed  to  enhance  the  spatial 
resolution  of  the  Landsat  Thematic  Mapper  thermal  band  from  the  remaining  bands.  The 
method  described  in  [2]  and  [7]  models  the  relationship  between  the  high-resolution  image 
and  the  image  to  be  estimated.  Regression  techniques  are  used  at  the  available  lower 
resolution  to  obtain  the  model  parameters. 

Another  approach  to  hyperspectral  resolution  enhancement  is  based  on  spectral  mixture 
analysis.  A  method  based  on  the  linear  mixing  model  has  been  investigated  that  employs 
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Figure  2.1:  Block  diagram  illustrating  the  flow  of  information  in  the  MAP  estimation  frame¬ 
work. 

constrained  nonlinear  optimization  techniques  to  obtain  high  resolution  endmember  frac¬ 
tions  [12-14],  Even  with  the  physical  constraints  that  assure  radiometric  validity  of  the 
derived  sub-pixel  spectra,  the  optimization  is  underdetermined  and  not  always  able  to  ar¬ 
rive  at  a  sharpened  solution.  An  alternative  approach  appends  the  high  resolution  image 
to  the  hyperspectral  data  and  computes  a  mixture  model  based  on  the  joint  data  set  [15]. 
A  high  resolution  hyperspectral  image,  however,  is  not  explicitly  estimated. 

Our  approach  is  based  on  a  novel  maximum  a  posteriori  (MAP)  estimation  frame¬ 
work  [17]  for  enhancing  the  spatial  resolution  of  an  image  using  co-registered  high  spatial- 
resolution  imagery  from  an  auxiliary  sensor.  The  estimation  framework  developed  allows 
for  any  number  of  spectral  bands  in  the  primary  and  auxiliary  sensor.  The  technique  is 
suitable  for  applications  where  some  correlation,  either  localized  or  global,  exists  between 
the  auxiliary  image  and  the  image  being  enhanced.  A  spatially  varying  statistical  model  is 
used  to  help  exploit  localized  correlation  between  the  primary  and  auxiliary  image.  Another 
important  aspect  of  the  algorithm  is  that  it  allows  for  the  use  of  an  accurate  observation 
model  relating  the  “true”  scene  with  the  low-resolutions  observations.  This  means  that 
a  potentially  wavelength-dependent  spatially- varying  system  point  spread  function  (PSF) 
can  be  incorporated  into  the  estimator.  Figure  2.1  illustrates  the  flow  of  information  in  the 
estimation  framework. 

Figures  2.2,  2.3,  2.4,  and  2.5  introduce  some  terminology  and  notation  used  throughout 
the  report.  Figure  2.2  shows  a  low  resolution  hypercube  consisting  of  M  =  mvrrih  spatial 
coordinates  and  P  spectral  coordinates.  When  represented  as  a  one  dimensional  vector,  a 
low  resolution  hypercube  is  denoted  by: 

y  =  [yi,  1,  •••)  2/p, i)  1/1,2,  •••;  Up, 2 -i  •••Wi.M,  -,Vp,m]T-  (2.1) 

Similarly,  Figure  2.3  shows  a  multispectral  data  cube  with  N  =  nvrih  spatial  coordinates 
(A  >  M)  and  v  spectral  coordinates  (v  <  P).  For  a  panchromatic  image,  v  —  1.  The  one 
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dimensional  vector  representing  a  multispectral  data  cube  is  denoted  by: 

X  [^1,1?  •••)  *£i/,l)  ^1,2,  •••)  *£z/,2 5  •••?  ^1,AT)  •••)  •  (2.2) 

High  resolution  hypercubes  are  denoted  by  2:  or  z.  As  depicted  in  Figure  2.4,  they  have 
N  spatial  coordinates  and  P  spectral  coordinates.  A  high  resolution  hypercube  in  one 
dimension  is  denoted  by: 

Z  =  [zi,i,  •••,  Zpt  1,  2:1,2, Zpt2, zi,w,  ^p,iv]T-  (2.3) 


Other  terminology  used  throughout  is  shown  in  Figure  2.5.  Note  that  high  and  low 
resolution  hyper-pixels  are  one  dimension  vectors  of  length  P  that  vary  in  the  spectral 
dimension  and  correspond  to  a  given  spatial  coordinate.  A  super-pixel  is  a  collection  of 
adjacent  high  resolution  hyper-pixels  that  exactly  fit  into  a  single  corresponding  low  reso¬ 
lution  hyper-pixel.  We  require  that  each  super-pixel  correspond  to  a  unique  low  resolution 
hyper-pixel  which  is  equivalent  to  —  and  —  being  integers. 

TTlv 

In  Equations  (2.1)  and  (2.3),  hyper-pixels  are  listed  in  column  order  and  similarly  for 
Equation  (2.2).  For  example,  the  first  nv  hyper-pixels  listed  in  Equation  (2.3)  represent  the 
first  column  of  hyper-pixels  shown  in  Figure  2.4. 


M  low  resolution  hyper-pixels 


Figure  2.2:  Low  Resolution  Hypercube  y 

Throughout  various  portions  of  the  report,  we  make  references  to  Mercury  computers, 
on  which  certain  CPU  intensive  parts  of  the  MAP  algorithm  were  implemented.  The  Mer¬ 
cury  compute  system  consists  of  multiple  compute  nodes  interconnected  by  a  high  speed 
interconnect  fabric  called  the  Raceway.  Each  compute  node  consists  of  a  high  speed  floating 
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N  high  resolution  multispectral  pixels 

Figure  2.3:  Multispectral  Data  Cube  x 

point  processor  with  an  attached  vector  processor  unit  (VPU).  Each  node  also  contains  local 
memory  and  an  interconnection  to  the  Raceway. 

This  report  contains  8  chapters  including  this  introduction.  Chapter  3  discusses  the 
delivered  software,  including  highlighted  features,  installation  procedure,  and  software  use. 
Chapter  4  presents  carefully  derived  explicit  formulas  directly  implemented  in  the  software 
provided.  It  also  answers  some  practical  questions  that  arose  during  the  initial  phases  of  im¬ 
plementation.  Chapter  5  analyzes  the  performance  of  the  algorithm  and  Chapter  6  presents 
additional  analysis  of  the  MAP  algorithm  with  emphasis  on  Mercury  timing.  Chapter  7 
discusses  Stochastic  Mixing  Models  and  provides  information  related  to  interfacing  them 
with  our  MAP  software.  Chapter  8  discusses  mis-registration  between  input  data  cubes  and 
its  effect  on  estimation  performance.  Appendices  A  and  B  detail  a  mathematical  analysis 
of  alternate  forms  of  the  MAP  algorithm,  and  Appendix  C  gives  a  atmospheric  analysis 
preformed  by  the  Nashua  office  of  ATK  Mission  Research. 
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N  high  resolution  hyper-pixels 


Figure  2.4:  High  Resolution  Hypercube  z 


High  resolution  hyper-pixel  zn 


Low  resolution  hyper-pixel  ym 


Super-pixel  zm 


Multispectral  pixel  xn 


Figure  2.5:  Notation  and  Terminology 
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Chapter  3 


Hyper  spectral  Image  Enhancement 
Software  Installation  and  Use 

A  main  concentration  of  the  Phase  II  effort  was  to  deliver  a  working  implementation  of 
the  MAP  algorithm  on  the  Mercury  Platform.  Since  one  of  the  goals  was  for  the  delivered 
software  to  be  easily  accessible  as  a  research  tool,  we  have  chosen  to  write  all  but  the 
most  CPU  intensive  portions  in  MATLAB.  Secondly,  due  to  the  potentially  long  run  times 
that  might  occur,  we  offload  the  CPU  intensive  parts  of  the  implementation  to  a  network 
of  compute  nodes  running  in  parallel.  Thirdly,  because  there  are  a  number  of  variations 
of  the  MAP  code,  we  incorporate  them  all  into  a  single  code  that  could  be  called  as  a 
function  from  MATLAB.  Lastly,  due  to  the  complexity  of  inputs  and  calling  arguments  and 
the  relationships  between  them,  we  developed  a  GUI  interface  that  will  reduce  the  time  it 
might  take  for  a  new  user  to  begin  using  the  software. 

We  want  to  emphasize  that  we  have  provided  a  powerful  tool  that  is  a  blend  of  code 
written  in  MATLAB  (for  ease  of  use  and  portability  across  computing  architectures)  and 
C  (for  speed  and  computing  power  on  a  Mercury  system).  The  use  of  MATLAB  as  a  user 
interface  drastically  simplifies  the  use  of  this  algorithm  and  relieves  the  user  from  having 
to  know  details  concerning  the  use  of  Mercury  systems.  All  the  user  has  to  know  is  how  to 
run  MATLAB  and  they  can  tap  into  the  power  of  the  Mercury  system. 

The  intensive  portions  of  the  MAP  Algorithm  are  implemented  on  a  Mercury  computer 
system.  However,  the  MATLAB  scripts  and  GUI  supplied  should  run  on  any  computing 
hardware  with  MATLAB  installed.1  The  user  only  need  set  the  number  of  available  Mercury 
CE’s  (Computational  Elements)  equal  to  0.  In  the  case  of  a  Mercury  computer  system  whose 
host  computer  has  MATLAB  installed,  any  number  of  available  CE’s  may  be  used.  Whether 
it’s  through  the  supplied  MAP_Algorithm.m  function  or  the  GUI  interface,  the  full  use  of  all 
nodes  on  a  Mercury  computer  has  been  made  simple  for  the  user. 

We  have  tested  the  software  on  four  Mercury  computer  systems  (two  with  a  Windows 
1The  GUI  requires  MATLAB  version  6.5  or  higher. 
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2000  PC  host  and  two  with  a  Unix  Sun  host)  and  as  a  result  are  confident  that  our  parallel 
implementation  is  portable  to  any  Mercury  system  with  MATLAB  installed  on  the  host 
computer.2  One  such  Mercury  system  (the  AdapDev  1280)  is  scheduled  for  delivery  at 
the  end  of  the  contract  and  will  come  with  all  required  software  pre-installed  (excluding 
MATLAB). 

As  discussed  below,  we  provide  a  single  MATLAB  function  MAP_Algorithm .  m  that  may 
be  directly  called  by  the  user.  Since  the  interface  is  through  MATLAB,  the  analyst  need 
not  be  concerned  with  the  details  of  parallel  programming,  but  can  still  benefit  from  the 
fast  response  times  resulting  from  it.  Other  advantages  of  using  MATLAB  include  plotting 
capabilities  and  a  familiar  computing  environment. 


3.1  Major  Functions,  Inputs,  and  Outputs 


Estimates  statistical  parameters 
and  loops  over  superpixels  serially. 


Estimates  statistical  parameters 
and  calls  Mercury  C  code 
parallelized  over  superpixels. 


Figure  3.1:  Major  Functions  of  the  MAP  Algorithm 

Figure  3.1  depicts  the  major  functions  of  the  MAP  Algorithm.  We  see  from  the  figure 
that  the  main  function  provided  is  MAP_Algorithm,  which  may  be  called  directly  by  the 
user  through  use  of  a  user  supplied  wrapper  function.  We  have  supplied  several  example 
wrapper  scripts,  including  an  easy  to  use  GUI.  The  two  other  main  functions  depicted 
compute  estimates  of  statistical  parameters  as  well  as  the  MAP  estimate.  The  “host” 
function  is  written  solely  in  MATLAB  and  is  portable  to  any  hardware  platform  running 
MATLAB.  The  host  function  makes  a  call  to  another  function  that  computes  the  MAP 
estimate  by  looping  over  super-pixels.  The  “mere”  function  runs  on  a  Mercury  platform.  It 
computes  the  statistical  estimates  in  MATLAB  and  the  MAP  estimate  through  a  parallel 
computation  over  super-pixels. 

Figure  3.2  depicts  the  major  inputs  and  outputs  of  MAP_Algorithm.  Perhaps  the  most 
important  inputs  are  a  low  resolution  hypercube  y  and  either  a  panchromatic  or  multispec- 

21  GB  or  more  of  physical  memory  on  the  host  computer  is  recommended. 
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Figure  3.2:  Major  Inputs  and  Outputs  of  the  MAP  Algorithm 


tral  image  array  x.  These  may  come  in  the  form  of  collected  data  from  a  two  sensor  platform 
or  as  degraded  images  from  a  known  high  resolution  hypercube  z.  The  main  output  is  the 
MAP  estimate  z  of  a  (possibly  hypothetical)  hypercube  z.  Some  of  the  other  inputs  of 
Figure  3.2  are  briefly  discussed  below. 

The  point  spread  function  is  linear  and  maps  a  hypothetical  vector  z  to  a  vector  y 
according  to  the  statistical  model  y  =  W z  +  n.  We  take  IF  to  be  a  weighted  averaging 
process  on  z.  The  entries  of  each  super-pixel  are  spatially  averaged  one  band  at  a  time 
in  order  to  obtain  the  corresponding  low  resolution  hyper-pixel.  Assuming  super-pixels 
do  not  overlap  spatially  yields  a  sparse  matrix  IF  with  orthogonal  rows  and  columns.  It 
may  be  partitioned  so  that  each  block  is  a  (possibly  0)  multiple  of  an  identity  matrix  of 
appropriate  size.  This  special  form  of  IF  facilitates  PCA  subspace  analysis,  because  the 
model  y  =  IFz  +  n  is  preserved. 

The  spectral  weight  matrix  is  optional  and  corresponds  to  a  similar  linear  model  x  = 
Sz  +  rj.  Intuitively,  this  model  should  be  assumed  when  frequencies  corresponding  to  x 
and  y  overlap  significantly  and  should  not  be  assumed  when  they  do  not.  Optionally,  the 
matrix  S  may  be  supplied  by  the  user  or  else  computed  from  x  and  y  using  a  least  squares 
technique. 
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We  have  provided  a  total  of  four  options  for  computing  statistical  parameters,  three  of 
which  are  self  contained,  and  one  of  which  requires  an  external  Stochastic  Mixing  Model 
computation.  The  statistical  parameter  estimation  plays  a  significant  role  in  the  determining 
the  quality  of  the  resulting  map  estimate  z.  Consequently,  we  have  anticipated  the  desire  to 
interface  new  parameter  estimation  algorithms  and  have  designed  the  software  accordingly. 

There  are  three  PCA  options  available.  The  first  option  is  to  not  use  PCA  at  all, 
but  rather  compute  solely  in  the  original  “spectral  components”.  The  second  option  is  to 
compute  the  MAP  estimate  in  a  PCA  vector  space  and  return  the  MAP  estimate  z  in  the 
same  space.  The  third  option  is  a  combination  of  the  first  two,  namely  to  compute  the 
MAP  estimate  in  PCA  space,  but  return  the  MAP  estimate  as  spectral  components.  Any 
number  of  PCA  components  may  be  chosen. 

Another  option  allows  the  application  of  a  hybrid  MAP-Interpolation  algorithm.  After 
specifying  the  number  of  components  to  be  used  for  MAP  estimation,  the  remaining  com¬ 
ponents  specified  are  generated  using  a  simple  interpolation  method.  Although  the  hybrid 
method  works  in  both  PCA  space  and  spectral  space,  it  is  probably  most  useful  in  conjunc¬ 
tion  with  PCA.  Generally,  there  is  little  to  gain  when  using  the  MAP  estimation  technique 
with  lower  ordered  PCA  components.  When  using  PCA,  a  hybrid  approach  typically  pro¬ 
duces  answers  that  compare  favorably  with  a  non-hybrid  approach,  but  runs  faster  than  the 
non-hybrid  approach. 


3.2  Software  Highlights 

As  shown  in  Figure  3.3,  MAP_Algorithm.m  logically  breaks  up  into  2  major  components 
depending  on  whether  the  underlying  broadband  frequencies  overlap  or  are  disjoint.  In  the 
former  case,  a  linear  model  x  =  Sz  +  77  between  the  broadband  (or  multispectral)  image 
and  the  unknown  high  resolution  hyperspectral  image  is  assumed;  in  the  later  it  is  not.  As 
detailed  in  Chapter  4  and  Appendices  A  and  B,  all  three  forms  of  the  MAP  algorithm  have 
been  analyzed  mathematically.  Under  the  assumption  of  a  linear  model  x  =  Sz  +  rj,  we 
have  proved  that  all  three  forms  of  our  MAP  estimator  are  identical  (Section  4.5).  If  no 
such  linear  model  is  assumed,  only  the  Form  1  MAP  estimator  is  applicable.  Therefore, 
the  Form  1  estimator  alone  is  sufficient  for  implementation,  although  it  requires  different 
implementations  depending  on  whether  the  linear  model  is  assumed  or  not. 

Each  implementation  above  breaks  down  further  into  a  General  Form  1  implementation 
and  a  Simplified  Form  1  implementation.  The  latter  assumes  a  more  restricted  structure  of 
the  associated  covariance  matrix  and  will  run  faster.  Thus,  there  are  4  major  components  of 
the  Form  1  MAP  estimator  that  have  been  implemented,  and  all  four  have  been  integrated 
into  the  single  user  interface  MAP_Algorithm.m. 

Referring  to  Figure  3.3,  the  fastest  run  times  are  realized  when  the  Simplified  Form  1 
Implementation  is  exercised.  This  occurs  when  a  simplified  zero- nonzero  structure  of  the 
covariance  matrix  is  assumed,  and  the  noise  parameter  an  is  zero.  The  simplified  structure 
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Unconditional  Covariance 
Estimate  Required. 


Conditional  Covariance 
Estimate  Required. 


Figure  3.3:  Logical  Break  up  of  the  MAP  Algorithm 


of  the  covariance  matrix  is  referred  to  as  the  “special  case  assumption” ,  and  the  requirement 
that  an  =  0  is  referred  to  as  the  “noise  free  condition” .  If  both  the  special  case  assumption 
and  the  noise  free  condition  hold,  a  matrix  inverse  is  avoided  for  each  super-pixel  resulting 
in  a  fast  MAP  estimation.  Assuming  a  fast  method  of  estimating  covariance  matrices  is 
available,  a  real  time  Mercury  implementation  is  feasible  based  on,  for  example,  Equation 
(4.70). 

The  Simplified  Form  1  implementation  is  automatically  detected  by  MAP_Algorithm.m. 
All  the  user  need  do  is  set  var_y=0  (aft  is  the  input  parameter  var_y  of  MAP_Algorithm .  m) 
and  specify  a  covariance  estimation  technique  consistent  with  the  special  case  assumption. 
Speed  improvements  are  automatically  enabled  when  the  noise  free  condition  and  special 
case  assumptions  are  encountered. 

In  general,  the  code  structure  assumes  all  covariance  matrices  are  block  diagonal,  where 
all  blocks  have  the  same  size,  but  may  otherwise  be  different  from  one  another.  However, 
the  special  case  assumption  requires  additional  structure.  While  the  blocks  need  not  be 
identical,  they  are  required  to  be  “identical  within  super-pixels”.  Details  are  presented  in 
Section  4.11  of  Chapter  4. 

As  seen  in  Figure  3.3,  the  covariance  estimate  required  is  unconditional  when  x  =  Sz  +  rj 
is  assumed,  and  is  conditional  otherwise.  In  either  case,  a  conditional  covariance  estimate 
is  ultimately  used  as  part  of  the  Form  1  MAP  Estimate.  If  the  linear  model  x  =  Sz  +  77  is 
assumed,  we  convert  it  to  the  conditional  version  using  Equation  (4.27).  This  is  desirable 
since  an  unconditional  covariance  is  easier  to  specify  than  a  conditional  covariance. 

There  are  presently  two  ways  to  specify  a  conditional  covariance  matrix  when  x  =  Sz+r] 
is  not  assumed.  Both  are  based  on  VQ  partitioning.  In  one  case  (ccov_meth=l),  VQ  clas- 
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sification  is  based  on  the  input  hypercube  y  and  in  the  other  case  (ccov_meth=2)  VQ  clas¬ 
sification  is  based  on  an  estimate  of  the  conditional  mean.  The  trade-off  between  these  two 
methods  is  speed  versus  image  enhancement  quality.  In  the  first  method  (ccov_meth=l), 
the  special  case  assumption  holds  so  that  no  matrices  are  inverted,  and  the  MAP  algorithm 
runs  quickly. 

There  is  presently  one  self  contained  algorithm  for  estimating  the  unconditional  covari¬ 
ance  matrix  when  x  =  Sz  +  r 7  is  assumed.  This  is  a  simple  estimation  technique  where  all 
blocks  of  the  block  diagonal  covariance  matrix  are  identical,  so  that  the  special  case  assump¬ 
tion  is  met  in  this  case  as  well.  The  Stochastic  Mixing  Model  (SMM)  interface  discussed 
later,  also  falls  under  the  assumption  x  =  Sz  +  77,  and  it  may  be  thought  of  as  a  second 
method  for  estimating  the  unconditional  covariance  matrix.  However,  this  method  is  not 
self-contained  since  it  requires  input  from  an  external  SMM  algorithm. 

We  have  established  a  software  framework  for  incorporating  new  covariance  estimation 
techniques.  A  MATLAB  programmer  can  incorporate  the  new  estimation  technique  by 
adding  a  new  function  call  to  the  existing  MATLAB  function  param_est  .m.  Any  new 
technique  that  adheres  to  the  general  input  and  output  data  structures  established  in 
param_est  .m  can  be  incorporated  in  a  straightforward  manner.  We  believe  that  these  data 
structures  are  sufficiently  general  for  the  easy  incorporation  of  new  covariance  estimation 
techniques  as  they  become  available. 

The  code  that  runs  on  a  Mercury  system  makes  extensive  use  of  the  Scientific  Algo¬ 
rithm  Library  (SAL),  which  is  a  C  library  provided  by  Mercury.  SAL  is  a  floating  point 
library  specifically  optimized  for  AltiVec  equipped  PowerPC  processors.  It  has  proven  to 
significantly  speed  up  codes  residing  on  a  single  node.  However,  the  functions  available  are 
somewhat  limited  as  compared  to  a  more  mature  library  such  as  LAPACK. 

We  also  make  extensive  use  of  another  C  library  provided  by  Mercury  called  the  Parallel 
Acceleration  System  (PAS).  The  PAS  library  is  a  means  by  which  the  available  nodes  pass 
data  to  one  another.  An  advantage  of  using  PAS  as  opposed  to  another  communication 
library  (e.g.  DX)  is  its  ability  to  scale.  For  example,  we  developed  our  C  code  on  the  8  node 
Mercury  AdapDev  1280,  and  then  we  successfully  ported  it  to  both  the  96  node  (Titan)  and 
64  node  (Hawk)  computer  systems  at  Wright- Patterson  Air  Force  Base  (WPAFB).  Due  to 
the  scalability  of  PAS,  the  port  was  accomplished  with  minimal  changes.  A  disadvantage 
of  PAS  is  that  its  portability  is  limited  to  Mercury  machines,  while  DX  is  standard. 

Through  use  of  the  PAS  communication  library,  we  have  developed  a  scalable  framework 
for  projects  implemented  on  Mercury  computer  systems.  The  intention  of  this  framework  is 
to  allow  programmers  with  little  or  no  knowledge  specific  knowledge  of  Mercury  computers 
to  develop  applications  relatively  quickly.  The  framework  also  facilitates  swapping  alterna¬ 
tive  functions  with  different  argument  lists  as  separate  building  blocks.  As  new  functions  are 
developed,  swapping  them  with  existing  functions  is  accomplished  with  minimum  impact 
to  the  PAS  framework  that  surrounds  them. 

Other  highlights  include  the  many  options  available  within  MAP_Algorithm.m,  some  of 
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which  are  listed  below: 


•  Work  in  either  the  spectral  domain  or  PCA  domain. 

•  Use  MAP  estimation  on  leading  components  and  interpolation  on  the  rest. 

•  Stochastic  Mixing  Model  interface. 

•  Flexibility  in  specifying  PSF  or  spectral  response  functions. 

•  Work  with  either  panchromatic  or  multispectral  input  data. 

These  options  are  discussed  further  in  the  documentation  below. 


3.3  Installation 

All  required  software  has  been  placed  in  a  single  file.  When  installed  (as  directed  below), 
the  result  is  a  subdirectory  named  mapalgorithm.  Once  you  have  generated  this  directory, 
all  of  the  required  files  will  be  available  to  MATLAB.  Windows  based  machines  require  the 
Winzip  utility  (a  scaled  down  version  of  Winzip  comes  with  Windows  XP),  and  Linux/Unix 
machines  require  the  tar  utility. 

For  Windows  based  machines,  unzip  the  compressed  file  mapalgorithm.zip  to  c:\ 
thereby  generating  the  directory  c :  \mapalgorithm.  For  Linux/Unix  based  machines,  untar 
the  mapalgorithm. tar  file  by  issuing  the  command  tar  xvf  mapalgorithm. tar. 

In  MATLAB,  change  the  working  directory  to  mapalgorithm  and  begin  running  the 
code  using  either  the  graphical  user  interface  or  the  provided  MATLAB  script  wrappers. 
For  Windows  based  machines,  the  byte  format  option  (specified  within  the  GUI  for  example) 
should  be  set  as  little  endian  (LE),  and  for  Linux/Unix  based  machines  it  should  be  set  as 
big  endian  (BE). 


3.4  Use  of  the  GUI 

Figure  3.4  displays  a  GUI  wrapper  of  MAP_Algorithm.m  named  HyperViewer  .m.  It  has  the 
same  functionality  as  the  MATLAB  script  MAP_Wrapper .  m  described  in  Section  3.5,  however, 
it  is  more  intuitive  and  simpler  to  use.  Certain  options  automatically  become  inaccessible 
as  appropriate  when  other  options  are  selected,  relieving  the  user  of  remembering  how  the 
various  options  relate  to  one  another. 

An  underlying  assumption  is  that  a  known  z  is  available  to  aid  us  in  determining  the 
performance  of  the  MAP  Algorithm.  This  known  z  will  be  used  to  generate  y  and  x.  Using 
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Figure  3.4:  HyperViewer  GUI 


x  and  y  as  input,  the  MAP  Algorithm  is  used  to  estimate  a  high  resolution  hypercube  z. 
Then,  z  and  z  are  compared  in  order  assess  the  effectiveness  of  MAP_Algorithm.m. 

The  starting  point  is  a  MATLAB  “.mat”  file  containing  a  high  resolution  hypercube, 
which  may  be  reduced  either  spatially  or  spectrally  in  order  to  form  z.  This  is  useful  if  a 
large  data  set  is  stored,  but  only  a  subset  of  it  is  desired  for  analysis.  The  “.mat”  file  must 
contain  the  “full  scale”  hypercube  as  a  three  dimensional  variable  with  vertical  dimension 
first,  horizontal  dimension  next,  and  spectral  dimension  last.  An  example  is  provided  in  the 
file  aviris.mat. 

After  loading  the  “.mat”  hie  with  the  MATLAB  Load  command,  HyperViewer  GUI  is 
started  by  typing  HyperViewer  (Z)  at  the  MATLAB  prompt  where  Z  is  the  name  of  the  full 
scale  hypercube  stored.  The  variable  name  supplied  may  be  anything,  but  it  is  identified  in 
HyperViewer  as  hyp.  In  Figure  3.4  we  see  the  full  scale  hypercube  being  displayed,  indicated 
by  the  drop  down  window  in  the  upper  left  hand  corner  showing  the  word  “hyp” . 

Once  the  MAP  Algorithm  is  run  by  clicking  the  Run  button,  the  scaled  down  hypercube 
z,  the  degraded  images  x  and  y,  and  the  MAP  estimate  z  all  become  available  for  display 
in  the  same  fashion  as  hyp.  The  low  resolution  hypercube  y  is  returned  as  two  possibly 
different  hypercubes,  one  (yin)  being  input  to  Map_Algorithm.m  and  the  other  (yOut)  being 
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output  from  Map_Algorithm.m.  They  will  be  different  only  if  PCA  analysis  is  specified  by 
taking  pcaMode=l. 

The  image  corresponding  to  any  band  may  be  selected  for  viewing  simply  by  selecting  a 
band  appearing  in  the  scroll  bar  immediately  to  the  right  of  the  image.  If  RGB  is  checked, 
then  three  bands  must  be  selected  in  order  to  produce  an  RGB  image.  This  is  done  by 
selecting  the  band  first  within  the  same  scroll  bar,  and  then  clicking  one  of  the  R,  G,  or  B 
buttons,  to  assign  the  band  to  Red,  Green,  or  Blue  respectively.  Alternately,  these  may  be 
entered  manually.  Two  dimensional  plots  of  band  number  versus  intensity  may  be  plotted 
by  checking  “zPlot” ,  and  then  clicking  on  the  desired  pixel. 

The  spatial  dimensions  of  the  initial  scaled  down  hypercube  z  are  specified  as  nv  and  rp, 
for  vertical  and  horizontal  respectively.  These  dimensions  must  be  less  than  or  equal  to  the 
dimensions  of  the  full  scale  hypercube  hyp.  The  spectral  dimension  of  z  (denoted  by  P )  is 
determined  by  specifying  Start,  Increment,  and  Stop  bands  of  hyp.  For  the  data  specified 
in  Figure  3.4,  the  bands  of  z  will  be  1,  6, 11, ... ,  221  so  that  P  =  45.  Thus,  z  is  simply  a 
restricted  version  of  the  full  scale  hypercube  hyp. 

The  creation  of  y  and  x  depend  on  the  dimension  parameters  rnv ,  rrp, ,  and  v  (displayed 
as  mv,  mh,  and  nu  resp.)  as  well  as  the  point  spread  function  PSF  and  the  spectral  response 
matrix  s.  The  high  resolution  nv  x  nn  x  P  hypercube  z  is  degraded  spatially  using  the  PSF 
to  obtain  the  mv  x  m/,  x  P  hypercube  y.  Similarly,  z  is  degraded  spectrally  using  s  in  order 
to  obtain  the  nv  x  rih  x  u  multispectral  data  cube  x.  The  user  has  the  flexibility  to  define 
his  own  PSF  or  spectral  response  matrix  s  by  editing  a  script.  Scripts  yielding  the  default 
PSF  and  matrix  s  are  provided  as  examples. 

The  PSF  used  to  determine  y  from  z  should  be  thought  of  as  a  two-dimensional  Finite 
Impulse  Response  (FIR)  filter.  Such  a  filter  may  be  represented  as  an  Lv  x  L/t  matrix  where 
Lv  =  nv/mv  and  =  nh/rrih .  The  spectral  response  function  s  is  represented  as  a  P  x  u 
matrix.  The  kth  column  of  s  consists  of  the  spectral  weights  used  to  determine  the  kth 
multispectral  value  of  x  via  weighted  averaging  of  high  resolution  hyper-pixels.  The  only 
general  requirements  for  s  are  that  each  entry  be  nonnegative  and  each  column  sum  be  1. 
Both  of  these  functions  are  discussed  in  detail  in  Chapter  4. 

After  selecting  a  PCA  mode  option,  a  covariance  estimation  technique,  the  number  of 
VQ  partitions  (if  applicable),  and  a  designation  of  whether  the  linear  model  relating  x  and 
z  is  assumed,  the  user  may  run  the  MAP  algorithm  by  clicking  the  run  button.  After  the 
run  is  complete,  an  SNR  plot  indicating  the  effectiveness  of  MAP_Algorithm.m  is  obtained 
by  clicking  the  SNR  button,  or  a  time  line  plot  may  be  displayed  if  numCE  >  0.  Example 
SNR  plots  are  presented  in  Chapter  5  and  example  time  line  plots  are  presented  in  Chapter 
6. 


After  clicking  the  run  button,  all  inputs  and  outputs  of  MAP_Algorithm.m  are  returned 
to  the  MATLAB  base.  The  user  may  then  execute  any  additional  MATLAB  functions  or 
scripts  at  hand.  For  example,  the  user  may  wish  to  issue  the  MATLAB  SAVE  command  in 
order  to  save  all  inputs  and  outputs  of  MAP_Algorithm.m  to  a  “.mat”  file.  Since  y  can  be 
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changed  by  MAP_Algorithm .  m,  input  and  output  versions  of  y  reside  in  different  structures. 
Similarly,  input  and  output  versions  of  s  reside  in  different  structures.  Under  the  file  menu 
of  the  HyperViewer  GUI,  the  user  may  also  print  a  screen  shot  of  the  current  state  of  the 
GUI  to  help  document  what  went  into  generating  outputs  of  MAP_Algorithm.m.  The  screen 
shot  looks  best  if  printed  in  landscape  format. 

Documentation  for  many  of  the  options  appearing  in  Figure  3.4  may  be  retrieved  by 
resting  the  mouse  cursor  over  the  corresponding  variable.  In  addition,  GUI  variables  are 
a  subset  of  the  input  variables  of  MAP_Algorithm.m,  which  are  described  in  detail  within 
Section  3.5  as  well  as  the  code  comments  of  MAP_Algorithm.m. 

Also  appearing  in  Figure  3.4  is  an  option  to  supply  a  Stochastic  Mixing  Model  (SMM) 
file.  The  format  of  the  required  file  is  detailed  in  Section  3.5,  and  a  general  description  SMM 
is  presented  in  Chapter  7.  An  SMM  file  named  Eismann.mat  is  included  as  an  example. 

The  byte  format  option  little  endian  (LE)  or  big  endian  (BE),  is  dictated  by  whether 
the  host  computer  is  operating  under  Windows  or  Linux/Unix  respectively.  It  must  be  set 
by  the  user  as  appropriate  whenever  a  Mercury  system  is  being  used.  The  input  cfg_f  ile 
is  described  in  the  next  section  and  also  applies  solely  to  Mercury  systems.  It  may  be  set 
for  the  HyperViewer  GUI  by  editing  HyperViewer  .m. 


3.5  Direct  Use  of  the  MAP  Algorithm  Function 

In  addition  to  the  HyperViewer  GUI,  we  have  created  the  MATLAB  script  MAP_Wrapper  .m. 
Both  HyperViewer  and  MAP_Wrapper  .m  may  be  viewed  as  alternative  wrappers  of  the  un¬ 
derlying  function  MAP_Algorithm.m.  The  wrapper  MAP_Wrapper . m  should  be  helpful  as  a 
guideline  to  users  who  wish  to  write  their  own  wrapper  of  MAP_Algorithm.m.  It  is  set  up  to 
make  changing  inputs  relatively  easy,  and  is  more  transparent  from  a  programmers  perspec¬ 
tive  than  HyperViewer.  The  script  is  run  merely  by  typing  MAP_Wrapper  at  the  MATLAB 
prompt. 

The  underlying  assumption  made  by  MAP_Wrapper  .m  is  the  same  as  that  made  by 
HyperViewer,  namely  that  we  begin  with  a  full  scale  hypercube  hyp  and  optionally  re¬ 
duce  it  to  a  sub- hypercube  named  z.  The  inputs  y  and  x  are  created  by  degrading  z 
through  use  of  the  PSF  and  s  respectively.  The  function  MAP_Algorithm.m  is  then  run  to 
produce  the  MAP  estimate  z,  which  is  compared  with  z  in  order  to  assess  the  effectiveness 
of  MAP_Algorithm.m. 

After  calling  the  function  MAP_Algorithm.m,  MAP_Wrapper . m  automatically  plots  images 
for  bands  1,2,3,  the  average  over  all  bands  of  z,  and  an  SNR  plot  comparing  z  with  z.  If 
z  has  been  transformed  into  its  PCA  components  (i.e.  pca_mode==l),  transforming  z 
into  its  PCA  components  is  required  before  a  meaningful  SNR  plot  can  be  made.  This 
transformation  is  performed  automatically  by  MAP_Wrapper  .m  as  needed.  Whenever  one  or 
more  CEs  are  specified,  MAP_Wrapper .  m  finishes  by  plotting  a  time  line  similar  to  the  output 
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of  Mercury’s  TATLView  program.  This  allows  for  a  detailed  timing  analysis  such  as  that 
presented  in  Chapter  6. 

Other  wrappers  that  might  be  considered  could  involve  different  underlying  assumptions. 
For  example,  in  a  realistic  situation,  the  user  would  not  have  a  high  resolution  hypercube 
(hyp)  to  start  with.  The  assumption  then  would  be  that  y  and  x  are  known  a  priori  and 
would  not  be  created  by  degrading  a  known  z.  Assuming  x  and  the  hypothetical  z  are 
related  through  a  linear  model  x  =  Sz  +  rj,  the  matrix  S  (created  from  the  smaller  matrix 
s )  may  not  be  known.  In  this  case,  the  user  may  choose  to  have  s  created  by  supplying  it 
to  MAP_Algorithm.m  as  an  empty  matrix. 

Optionally  in  MAP_Algorithm.m,  a  hybrid  MAP/interpolation  may  be  specified  where 
the  first  Phyb  components  are  computed  using  the  MAP  estimate  and  the  last  P  —  Phyb 
are  interpolated.  For  simplicity,  this  option  is  not  exercised  in  either  HyperViewer  or 
MAP_Wrapper  .m,  however  it  has  been  thoroughly  tested  and  works  with  either  spectral  or 
PCA  components.  The  utility  of  the  hybrid  option  is  that  it  is  very  fast  compared  to  the 
general  MAP  algorithm  where  many  P  x  P  matrices  are  inverted.  Of  particular  interest  is 
the  hybrid  option  in  conjunction  with  PCA  analysis,  since  interpolation  on  higher  ordered 
components  gives  essentially  the  same  SNR  as  the  MAP  estimate. 

We  conclude  this  section  by  describing  all  input  and  output  arguments  of  MAP_Algorithm .  m. 
As  of  the  writing  of  this  document,  the  description  below  matches  the  code  comments  ex¬ 
actly. 


7. 

7.  Usage:  mapStructOut  =  MAP_Algorithm  (mapStructln) 

7. 

7.  The  I/O  arguments  to  map  algorithm  are  in  the  form  of  Structures  where 
7.  the  fields  of  the  structures  are  given  below  (i.e.  y  =  mapStructln. y)  . 

7. 


7,  Input 
7.  y: 

7. 

7o  x: 

7. 


7. 


7. 


7. 

7. 


(mapStructln) : 

Low  resolution  hypercube  (mv  x  mh  x  P) . 


Co-registered  and  aligned  high-resolution  multispectral  cube 
(nv  x  nh  x  nu) .  Here  aligned  means  that  x  and  y  span  the 
*exact*  same  field  of  view  and  nv/mv  and  nh/mh  are  integers. 
Such  alignment  is  crucial  to  performance.  For  a  panchromatic 
image ,  nu  =  1 . 


7«  hybrid:  1  to  apply  MAP/interpolation  hybrid,  0  for  MAP  alone. 

7. 


7.  P_hyb : 
7. 

7. 

7. 


Number  of  components  used  in  MAP  portion  of  hybrid.  The 
remaining  P  -  P_hyb  components  are  3d  interpolated. 
Requirements:  0  <=  P_hyb  <  P  if  hybrid  =  1, 

P_hyb  <=  Npca  if  hybrid  =  1  &  pca_mode  >  0. 
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7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


Set  P_hyb  =  0  for  pure  interpolation  (no  MAP  estimation) . 

pca_mode:  >  0  for  applying  MAP  algorithm  to  principal  components,  0 

for  applying  MAP  algorithm  to  original  spectral  components. 

Set  pca_mode=l  to  return  PCA  coordinates  and  set  pca_mode=2 
to  return  spectral  coordinates. 

Npca:  Number  of  principal  components  computed  if  pca_mode  >  0. 

Requirements:  1  <=  Npca  <=  P  if  pca_mode  >  0, 

P_hyb  <=  Npca  if  hybrid  =  1,  pca_mode  >  0. 

SMM_file:  A  MATLAB  ’.mat’  file  containing  high  res  abundance  maps, 
end-member  mean  vectors,  end-member  covariances,  and  a 
PCA  transformation  matrix  if  pca_mode  ==  1 .  If  there  is 
no  SMM  file,  set  SMM_file  =  ’none’.  An  SMM  file  is  expected 
to  contain  the  following  variables  with  the  exact  names  and 
sizes  as  specified  below  (also  see  MAP_SMM.m): 

a_smm  =  High  resolution  abundance  mapping  (nv  x  nh  x  Ne) 
where  Ne  is  the  number  of  end  members. 
m_eps  =  End  member  mean  vectors  (K  x  Ne)  where  K  =  P  if 
pca_mode  =  0,  and  K  =  Npca  otherwise. 

C_eps  =  End  member  covariances  (K  x  K  x  Ne) 

E_smm  =  PCA  transformation  matrix  (P  x  Npca) .  Take  E_smm 
as  an  empty  matrix  if  no  PCA  analysis  was 
performed,  in  which  case  K  =  P  and  a_smm,  m_eps, 
C_eps  are  all  ’spectral  domain’  quantities. 

Caution:  If  an  SMM  file  is  provided,  other  input  parameters  to 
this  function  must  be  provided  in  a  consistent  manner.  Of 
particular  concern  is  the  consistency  of  x,  y,  s,  and  nu.  In 
general,  whatever  (overlapping)  input  parameters  used  in  the 
creation  of  a_smm,  m_eps,  C_eps,  and  E_smm,  should  be  exactly 
the  same  parameters  provided  to  this  function. 

Jseed:  Initial  random  number  seed  for  VQ  partitioning. 

var_y:  Noise  variance  for  all  pixels  and  bands  in  y.  In  the  final 

report,  this  is  defined  as  the  variance  of  n  =  y  -  Wz.  It 
is  also  a  diagonal  load  of  an  inverted  matrix  appearing  in 
the  MAP  estimate.  Increasing  var_y  will  improve  the 
conditioning  of  the  matrix  being  inverted. 

var_x:  Noise  variance  for  all  pixels  and  bands  in  x.  In  the  final 
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report,  this  is  defined  as  the  variance  of  eta  =  x  -  Sz.  It 
is  also  a  diagonal  load  of  an  inverted  matrix  appearing  in 
the  equation  that  relates  conditional  parameter  estimates 
(which  are  required  for  our  Form  1  implementation)  with  the 
corresponding  unconditional  parameters.  Increasing  var_x 
will  improve  the  conditioning  of  the  matrix  being  inverted. 

psf :  FIR  point  spread  function  of  size  Lv  x  Lh  where  Lv  =  nv/mv, 

and  Lh  =  nh/mh. 

As  described  in  the  final  report,  specification  of  a  psf  matrix  of  size 
Lv  x  Lh  is  equivalent  to  taking  all  rows  of  the  M  x  L  matrix  omega_hat 
to  be  psf ( : ) ’  where  M  =  mv*mh  and  L  =  Lv*Lh.  As  presented  in  the 
report,  omega_hat  contains  the  non-zero  entries  of  the  M  x  N  matrix 
omega  and  the  MP  x  NP  matrix  W  is  given  in  terms  of  omega  by: 

W  =  Q*kron(omega,eye(P)) . 

xeqSz:  1  if  x  and  z  are  linearly  related  by  x  =  Sz  +  eta. 

0  otherwise. 

cov_meth:  Method  of  estimating  the  (unconditional)  covariance  matrix  Cz 
This  method  applies  only  when  xeqSz=l.  Presently  one  type: 
cov_meth  =  1:  all  PxP  blocks  are  identical. 

ccov_meth:  Method  of  estimating  the  conditional  covariance  matrix  Cz|x 
This  method  applies  only  when  xeqSz  =  0.  Presently  two  types: 
ccov_meth  =  1 :  VQ  partitioning,  classification  on  y. 
ccov_meth  =  2:  VQ  partitioning,  classification  on  Muz|x. 

I_meth:  Method  of  interpolation  for  obtaining  mean  images  from 

down  sampled  images.  Same  options  as  interp2  (e.g. ’linear’) . 

NVQ:  Number  of  partitions  in  VQ  (if  applicable). 

s:  Spectral  response  matrix  if  applicable  (i.e.  if  xeqSz  =  1). 

Specify  as  an  empty  matrix  in  order  to  compute  it  using 
a  constrained  least  squares  technique.  If  nonempty,  the 
size  of  s  should  be  P  x  nu  unless  both  hybrid=l  and 
pca_mode=0,  in  which  case  the  size  should  be  P_hyb  x  nu. 

Each  entry  of  s  must  be  nonnegative  and  each  column  must 
sum  to  unity. 

The  P  x  nu  spectral  response  matrix  s  (if  provided  at  all) 
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7. 

7. 

7. 

7«  num_CE: 
7. 

7. 

7. 

7. 


should  be  provided  in  the  spectral  domain.  If  PCA  analysis 
is  to  be  performed  (i.e.  pca_mode  >  0) ,  s  will  be 
automatically  transformed  into  PCA  space.  Therefore,  if 
pca_mode  >  0,  s  will  be  transformed  by  s  =  Ey’*s,  where  Ey 
is  the  P  x  Npca  matrix  described  below.  The  resulting  size 
of  s  is  Npca  x  nu. 

If  hybrid  =  1  and  pca_mode  >  0,  the  Npca  x  nu  matrix  s  is 
later  restricted  to  its  first  P_hyb  components.  Thus, 
regardless  of  the  value  of  pca_mode,  the  final  size  of 
s  is  P_hyb  x  nu  whenever  xeqSz  =  hybrid  =  1. 

Number  of  nodes  (CE’s)  to  use  on  a  Mercury  computer  system. 
Set  this  input  to  0  if  not  running  on  a  Mercury.  Otherwise, 
setting  this  to  0  will  run  MAP_Algorithm  on  whatever  host 
computer  is  being  used  as  interface  to  the  Mercury. 


7«  Input  required  if  running  on  a  Mercury  computer  system: 

7«  WorkerStack:  size  of  the  stack  on  each  worker  CE  (#  bytes) 

7«  WorkerHeap:  size  of  the  heap  on  each  worker  CE  (#  bytes) 

7o  ControlStack:  size  of  the  stack  on  the  controller  CE  (#  bytes) 

7o  ControlHeap:  size  of  the  heap  on  the  controller  CE  (#  bytes) 

7«  The  controller  heap  size  may  need  to  increase  when  the  controller 

7«  is  the  only  node  specified  for  execution.  It  may  be  decreased 

7«  proportionally  as  worker  nodes  are  specified,  because  each  node 

7o  will  need  proportionally  less  memory  to  hold  its  assigned 

7o  super-pixels.  Page  faults  and  resources  exhausted  errors  can 

7o  many  times  be  corrected  by  adjusting  the  controller  heap  size. 

7«  ceid_wkrs :  array  of  worker  nodes  to  run  on  the  Mercury  system 

7«  ceid_ctrl :  the  CE  id  of  the  control  node  on  the  Mercury  system 

7o  byte_format:  Byte  format  of  host:  1  =  little  endian  (e.g.  Intel  host) 

7o  2  =  big  endian  (e.g.  Sun  host) 

7o  cfg_file:  Name  and  location  of  Mercury  ’.cfg’  file  for  application 
7«  of  a  configmc  command  on  the  Mercury  host.  For  example, 

7«  cfg_file  =  ’  C:\MercurySoftware\AdapDev\AdapDev.cfg’ 

7o  will  issue  the  following  command  from  MAP_estimate_merc .m: 

7o  configmc  -cf  C:\MercurySoftware\AdapDev\AdapDev.cfg  init 

7o  The  intent  is  to  have  an  automatic  way  of  avoiding  "out  of 

7«  resources"  crashes  on  the  Mercury.  The  variable  cfg_file  is 

7«  ignored  if  not  running  on  a  Mercury  (designated  by  taking 

7«  numCE  =  0)  .  To  avoid  the  application  of  a  configmc  command 

7o  when  running  on  a  Mercury,  set  cfg_file  =  ’skip’. 

7. 
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7. 

7o  Output  (mapStructOut)  : 

7. 


7o  Pzhat : 

7. 

7. 

7. 

7.  y: 

7. 

7. 

7. 

7.  s: 

7. 

7. 

7. 

7. 

7o  zhat : 

7. 

7. 

7. 

7. 

7. 

7.  Cy : 

7. 

7. 

7. 

7.  Ey: 

7. 

7. 

7. 

7. 


Size  of  the  spectral  dimension  of  the  MAP  estimate.  It  is 
identical  to  the  input  P  (i.e.  third  dim  of  input  y)  unless 
pca_mode  =1.  If  pca_mode  =  1,  then  Pzhat  =  Npca. 

Low  resolution  hypercube  of  size  mv  x  mh  x  Pzhat.  If 
pca_mode  =  1,  y  is  transformed  into  PCA  space.  Otherwise,  it 
is  identical  to  the  input  y. 

Spectral  response  matrix  (if  applicable)  adjusted  according  to 
a  principal  component  transformation  if  pca_mode  =  1.  The 
output  size  is  Pzhat  x  nu  unless  hybrid  =  1 ,  in  which 
case  the  output  size  is  P_hyb  x  nu. 

The  MAP  estimate:  a  high  resolution  hypercube  of  size 
nv  x  nh  x  Pzhat.  If  pca_mode  =  1,  principal  components  are 
returned.  If  pca_mode  =  2,  processing  is  performed  in  PCA 
coordinates  but  the  resulting  MAP  estimate  is  returned  in  the 
original  spectral  coordinates. 

P  x  P  spectral  covariance  matrix  used  for  PCA  transformation. 
Returned  as  an  empty  matrix  if  pca_mode  =  0.  Here,  P  is  the 
size  of  the  spectral  dimension  of  the  input  y. 

P  x  Npca  matrix  of  eigenvectors  of  Cy.  The  jth  column  of  Ey 
is  the  eigenvector  corresponding  to  the  jth  eigenvalue  of  Cy 
where  the  Npca  eigenvalues  sorted  from  largest  to  smallest. 
Returned  as  an  empty  matrix  if  pca_mode  =  0. 


3.6  Other  Wrappers 

Also  provided  is  a  third  wrapper  of  MAP_Algorithm.m,  named  MAP_Wrapper_di agnostic  .m, 
whose  purpose  is  to  provide  the  user  with  diagnostic  information  about  the  magnitude  of 
differences  between  running  on  the  host  and  running  on  the  network  of  CEs.  An  important 
application  of  this  is  to  identify  malfunctioning  nodes.  We  have  used  the  diagnostic  wrapper 
to  identify  a  bad  node  (node  77  as  determined  by  the  default  configuration  file)  on  the  Titan 
96  node  Mercury  system  at  WPAFB.  The  cause  is  likely  to  be  faulty  local  memory  at  that 
node. 
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Code  run  on  the  network  of  CEs  is  written  in  single  precision  C.  Conversely,  code  run  on 
the  host  computer  is  written  in  double  precision  MATLAB.  Due  to  the  different  precisions 
involved,  there  will  be  differences  between  the  resulting  z  that  each  produces.  In  order  to 
quantify  these  differences,  MAP_Wrapper_diagnostic  .m  calls  the  MAP_Algorithm  twice  in 
succession:  first  on  the  host  computer,  and  the  next  on  the  network  of  CEs.  Afterward,  a 
plot  is  produced,  showing  the  magnitude  of  the  (mean)  differences  on  the  vertical  axis  and 
the  CE  identification  on  the  horizontal  axis. 

Editing  the  script  MAP_Wrapper_diagnostic  .m,  the  user  may  specify  any  desired  input 
parameters.  The  value  of  mapStructIn.num_ce  is  automatically  set  to  zero  for  the  first 
call  to  MAP_Algorithm.m,  and  automatically  set  to  the  user  specified  value  of  the  variable 
num_ce  for  the  second  call.  After  the  desired  inputs  are  supplied,  the  diagnostic  wrapper  is 
run  as  a  MATLAB  script  simply  be  typing  MAP_Wrapper_diagnostic  .m  at  the  MATLAB 
prompt. 

Figure  3.5  shows  an  example  of  running  MAP_Wrapper_diagnostic  .m  on  the  Mercury 
AdapDev  1280.  Each  marker  on  the  plot  represents  the  difference  in  z  corresponding  to  a 
particular  super-pixel.  The  CE  that  processed  the  super-pixel  lies  directly  underneath  the 
marker.  We  see  in  Figure  3.5  that  all  differences  between  the  two  runs  are  relatively  small, 
indicating  that  both  implementations  agree  with  one  another  and  that  all  nodes  tested  are 
functioning  properly.  An  unusually  large  difference  at  a  particular  CE  would  indicate  a 
problem  at  that  CE. 

A  fourth  wrapper  of  MAP_Algorithm . m,  named  plot_timing.m,  generates  timing  plots 
as  a  function  of  available  nodes.  In  this  wrapper,  several  runs  of  MAP_Algorithm.m  are 
made  is  succession,  where  the  number  of  nodes  used  for  a  given  run  varies  while  the  other 
input  variables  are  fixed.  These  plots  help  access  the  utility  of  the  Mercury  platform  as  it 
pertains  to  this  project.  Several  examples  are  presented  in  Chapter  6. 
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Comparson  of  Matlab  and  Mercury  generated  zhat 


Figure  3.5:  Diagnostics  on  the  AdapDev  1280 
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Chapter  4 

Mathematical  Derivations  and 
Analysis 


This  chapter  presents  the  mathematical  framework  relating  to  the  theory  and  computer 
implementation  of  a  Maximum  A  Posteriori  (MAP)  algorithm  for  image  enhancement.  A 
number  of  practical  questions  that  arose  during  the  initial  phases  of  implementation  are 
addressed  and  answered  here.  Some  of  the  topics  addressed  are  the  singular  or  non-singular 
status  of  certain  matrices  that  require  inversion,  the  implication  of  making  various  covari¬ 
ance  matrices  diagonal,  the  relationship  between  conditional  independence  and  diagonal 
covariance  matrices,  and  the  conditions  under  which  the  various  forms  of  the  MAP  algo¬ 
rithm  are  identical.  Careful  attention  has  been  paid  to  accurately  detail  the  mathematical 
equations  required  for  implementation  on  either  a  serial  or  parallel  computer  system. 

Almost  universally,  we  assume  the  same  notation  that  appears  in  the  article  [18].  In 
particular,  a  broadband  image  array  x  (which  is  generalized  to  a  multispectral  image  in 
Sections  4.14  and  4.13),  a  low  resolution  hypercube  y ,  and  a  high  resolution  hypercube  z 
have  entries  in  band-sequential  lexicographical  order: 

X  =  (x1,x2,...,xn)t, 

U  =  (2/1,1,  2/2,1,  •  •  • ,  Up, i,  2/1,2, 2/2,2,  •  •  • ,  2/p, 2,  •  •  • ,  2/i,m,  2/2, m,  •  •  • ,  2 /p,m)T, 

z  =  (zi,l,  22,1)  •  •  •  ,  zP,li  A, 2,  ^2,2,  •  •  •  ,  Zp, 2,  •  •  •  ,  2l,JV,  22,1V >  •  •  •  ,  2p,lv)T- 

Thus,  a  broadband  image  array  consists  of  N  pixels,  a  high  resolution  hypercube  con¬ 
sists  of  P  bands  with  N  pixels  in  each  band,  and  a  low  resolution  hypercube  has  P 

th 

bands  with  M  <  N  pixels  in  each  band.  The  m  low  resolution  hyper-pixel  is  given 
by  Vm  =  (2/1, rn,  2/2, rn,  •  •  • ,  2/P,m)T,  and  the  high  resolution  hyper-pixel  is  given  by  z)  = 
(zi,i,Z2,i,  ■  ■  ■  ,zPti)T.  A  super-pixel  is  a  collection  of  adjacent  high  resolution  hyper-pixels 
and  are  associated  with  a  given  low  resolution  hyper-pixel  as  described  in  Section  4.2. 

As  seen  in  the  above  ordering  for  y  and  z.  the  spectral  dimension  of  length  P  is  listed 
first.  We  denote  the  vertical  and  horizontal  spatial  dimensions  of  y  by  mv  and  m/,  and 
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those  of  z  by  nv  and  n/,  so  that  M  =  mvrrih  and  N  =  rq,n/,.  Referring  to  Figure  4.1  (where 
nv  =  4,  rih  =  6,  rnv  =  2,  and  m/(  =  3),  we  see  that  the  high  resolution  hyper-pixel  z%  lying  in 
row  k  and  column  j  has  index  i  —  (j  —  \)nv  +  k  for  1  <  k  <  nv  and  1  <  j  <  rif,.  Similarly, 
the  low  resolution  hyper-pixel  ym  lying  in  row  k  and  column  j  has  index  m  =  (j  —  l)rnv  +  k 
for  1  <  k  <  mv  and  1  <  j  <  rn/, .  Therefore,  hyper-pixels  y±,y2, ,  ijm  and  zlt  z2, . . . ,  zn 
are  listed  in  column  order.  The  entries  of  the  broadband  image  array  x  are  listed  in  column 
order  as  well. 

We  assume  that  z  is  a  multivariate  random  variable  with  a  prior  Gaussian  distribution 
and  further  assume  Bayesian  general  linear  models  (see  10.6  of  [17]): 

x  =  Sz  +  rj, 
y  =  Wz  +  n, 

where  S  denotes  an  N  x  NP  spectral  response  matrix1,  and  W  denotes  the  MP  x  NP 
matrix  form  of  a  point  spread  function.  By  definition  of  a  Bayesian  general  liner  model,  the 
vectors  r)  and  n  are  zero  mean  multivariate  Gaussian  random  variables  and  are  independent 
of  z.  It  follows  from  these  assumptions  that  x  and  z  are  jointly  Gaussian  as  are  y  and  z. 

The  model  x  =  Sz+y  is  only  assumed  optionally,  since  we  also  want  to  examine  the  case 
where  it  may  not  hold.  The  model  would  most  naturally  be  assumed  when  the  frequency 
band  associated  with  the  known  hypercube  y  contains  the  frequency  band  associated  with 
the  panchromatic  image  x.  The  model  would  most  likely  not  be  assumed  if  these  frequency 
bands  are  disjoint.  Much  of  this  chapter  applies  regardless  of  whether  or  not  the  model  is 
assumed.  The  way  to  decern  if  an  equation  depends  on  the  model  is  by  simply  observing 
whether  the  matrix  S  (or  its  related  quantity  s)  appears. 

Although  some  results  are  presented  in  terms  of  arbitrary  covariance  matrices  Cn  and 
Cn  of  rj  and  n,  we  will  typically  make  the  additional  assumption  that  they  are  multiples  of 
the  appropriate  identity  matrix: 


Cij  G7  In  i 

Cn  =  cr^IMP. 

There  are  three  equivalent  forms  of  MAP  estimators  address  here.  Each  form  is  rep¬ 
resented  in  terms  of  a  product  of  probability  density  functions  that  is  to  be  maximized 
over  all  possible  z.  The  argument  z  for  which  maximum  is  attained,  is  expressed  as  the 
solution  to  a  matrix  equation.  The  derivation  of  each  such  matrix  equation  from  its  PDF 
representation,  and  a  proof  of  their  equivalence,  is  presented  in  Section  4.5. 

1The  way  S  is  defined  later  in  Equation  (4.43),  x  —  rj  =  Sz  will  provide  x  —  rj  in  super-pixel  ordering 
(defined  later),  assuming  z  is  given  lexicographically. 
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4.1  Kronecker  Products 


Since  significant  use  is  made  of  the  direct  sum  and  the  Kronecker  tensor  product  [19],  we 
define  them  below  and  list  some  of  their  properties.  The  direct  sum  of  any  matrices  A  and 
B  is  given  by: 

A®B  =  diag  (A,B)  = 

Extending  this  in  the  obvious  way  to  a  direct  sum  of  M  matrices  we  have: 


A  0 
0  B 


'  M  \  -1  M 

m=l  /  m= 1 


where  Am  and  Bm  are  assumed  to  have  compatible  sizes  for  matrix  multiplication. 

For  any  matrices  A  and  B,  the  (right)  Kronecker  tensor  product  of  A  and  B  is  defined 
by: 


A®  B 


ai^B  . . 

CLlnB 

&2,1  B 

^2  ,nB 

^ .  1  B 

&m,nB 

where  ahJ  denotes  the  (i,j)  entry  of  the  m  x  n  matrix  A.  The  Kronecker  product  is  an 


intrinsic  MATLAB  function  “kron”  where: 

A®B  =  kron(A,  B). 

Several  useful  facts  concerning  Kronecker  tensor  products  are: 

(. A®B)®C  =  A®  (B  ®  C),  (4.1) 

A®(B  +  C)  =  A®  B  +  A®  C,  (4.2) 

(B  +  C)®A  =  B  ®  A  +  C  ®  A,  (4.3) 

(. A®B)t  =  At®Bt,  (4.4) 

(Ai  ®  -Bi)(A.2  0  B2)  =  A1A2  ®  B1B2 ,  (4-5) 


where  the  matrices  appearing  within  ordinary  matrix  operations  have  compatible  sizes.  The 
additional  property  ( A®v)B  =  AB  ®  v  where  v  is  a  column  or  row  vector,  may  be  derived 
from  (4.5)  by  setting  B  =  B  ®  1,  and  similarly  for  the  property  (v  ®  A)B  =  v  ®  AB. 
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4.2  Super-pixel  Ordering 


Figure  4.1:  Hypercube  partitioned  into  super-pixels 


Figure  4.1  depicts  a  high  resolution  hypercube  with  N  =  24,  M  =  6,  and  P  =  4 
partitioned  into  super-pixels.  The  column  ordering  of  the  high  resolution  hyper-pixels  is  z  = 
(zi;  z2; . . . ;  Z24),  where  the  semicolon  denotes  vertical  concatenation.  The  above  partitioning 
defines  a  new  ordering  given  by  listing  super-pixels  in  column  order.  Defining  q  to  be  the 
permutation 


/  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  \ 
y  1  2  5  6  3  4  7  8  9  10  13  14  11  12  15  16  17  18  21  22  19  20  23  24  )  ’ 


the  super-pixel  ordering  of  z  is  given  by  z  =  (zq(i);  zq^)]  •  •  • ;  zg( 24)).  Note  that  the  ordering 
of  y  is  unaffected  since  super-pixels  are  listed  in  column  order.  These  remarks  extend  to 
the  general  case  where  N  =  nhnv  and  M  =  mhmv  provided  that  77-  and  77-  are  integers. 
In  this  case,  each  super-pixel  is  a  ^x^xP  hypercube  and  the  hyper-pixels  of  adjacent 
super-pixels  do  not  overlap.  In  the  above  example  nv  =  4,  rif,  =  6,  rnv  =  2,  and  m/,  =  3.  The 
subscripts  h  and  v  are  used  to  represent  horizontal  and  vertical  directions  of  a  hypercube. 

The  NP  x  NP  permutation  matrix  corresponding  to  q  is  obtained  by  partitioning  the 

th 

NP  x  NP  identity  matrix  into  P  x  P  blocks  and  for  1  <  j  <  N,  placing  its  j111  block 
column  into  block  column  number  q(j).  Denoting  the  resulting  permutation  matrix  by  Q, 
the  high  resolution  hypercube  in  the  super-pixel  ordering  is  simply  z  =  Qz.  Since  Q  is  an 
orthogonal  matrix,  Q~L  =  QT  will  always  hold.  In  the  above  example,  we  have  in  addition 
Q  =  QT ,  however,  this  is  not  true  in  general  (the  case  (nv,rih,mv,mh)  =  (9,2,3, 1)  is  a 
counterexample).  For  programming  efficiency,  any  matrix  operation  involving  Q  may  be 
replaced  by  a  corresponding  operation  with  q. 
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4.3  Point  Spread  Functions 


The  matrix  W  appearing  in  the  Bayesian  general  linear  model  y  =  Wz  +  n  represents  the 
Point  Spread  Function  (PSF).  It  is  an  MP  x  NP  matrix  that  when  applied  to  z  performs 
blurring  and  down-sampling.  For  the  purposes  of  discussion  and  introduction,  we  assume 
for  the  moment  that  W  represents  an  un-weighted  averaging  transformation:  the  high 
resolution  hyper-pixels  belonging  to  a  super-pixel  are  averaged  to  obtain  the  corresponding 
low  resolution  hyper-pixel.  This  will  subsequently  be  generalized  to  a  weighted  average.  In 
terms  of  the  permuted  high  resolution  hyper-pixel  z,  the  linear  model  is  y  =  WQTz  +  n, 
or  equivalently  y  =  Wz  +  n  where  W  =  WQT . 

Since  both  W  and  W  are  MP  x  NP  matrices,  they  may  partitioned  into  P  x  P  blocks 
where  each  block  acts  on  a  high  resolution  hyper-pixel.  Since  W  represents  an  averaging  of 
super-pixels  and  there  are  N/M  hyper-pixels  per  super-pixel,  each  block  is  either  the  zero 
matrix  or  else  is  jflp.  The  nonzero  blocks  of  W  have  a  simpler  form  than  those  of  W.  For 
the  previous  example  of  IV  =  24,  M  —  6,  and  P  =  4,  the  matrix  W  has  the  form: 


W 


*  *  ** 


*  *  ** 


0 


*  *  ** 


0 


*  *  ** 


*  *  ** 


*  *  ** 


(4.6) 


where  each  *  represents 

This  simple  form  is  preserved  whenever  super-pixels  are  non-overlapping.  When  this 
occurs,  each  row  of  W  consists  of  ||  nonzero  P  x  P  blocks  j-Ip  placed  diagonally  as  in 
Equation  (4.6).  Now  let  L  —  ^  and  define  Wp  to  be  the  P  x  LP  matrix  consisting  of  L 
copies  of  j-Ip  concatenated  horizontally.  Then,  Wl  may  be  succinctly  denoted  by: 

WL  =  jJlxL®IP,  (4.7) 

where  Jixl  is  the  1  x  L  vector  of  all  ones  and  0  is  the  Kronecker  tensor  product.  Likewise, 
the  point  spread  function  W  is: 


M 

w  =  Im®Wl  =  Q)Wl.  (4.8) 

m= 1 

Thus,  W  consists  of  M  non  square  blocks  on  the  diagonal  (each  of  size  P  x  LP)  resulting 
in  a  matrix  of  size  MP  x  NP.  We  will  refer  to  this  example  as  the  un-weighted  averaging 
PSF. 

We  now  relax  the  assumption  that  W  has  the  form  (4.8)  but  still  insist  on: 
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1.  Super-pixels  do  not  overlap. 

2.  The  PSF  does  not  vary  spectrally. 


As  we  have  seen,  the  first  assumption  permits  the  use  of  a  permutation  q  to  describe  the 
super-pixel  ordering  of  high  resolution  hyper-pixels.  The  second  assumption  means  that  the 
PSF  is  fully  prescribed  by  an  arbitrary  M  x  N  matrix  u>.  If  we  fix  a  specific  spectral  value  A 
from  among  the  P  available,  and  consider  the  N  vector  zy  and  M  vector  y\  corresponding 
to  it,  then  y\  =  luz\  +  n\  where  n,\  is  the  associated  noise.  The  meaning  of  the  second 
assumption  is  that  u  does  not  depend  on  A. 

For  general  u>,  the  PSF  matrix  IF  is  given  by: 

IF  =  (jj  ®  Ip. 


In  particular,  from  Equations  (4.7)  and  (4.8), 


W  =  I 


M 


1 


Jl 


xL 


Im  <H>  JlxL 


I  Pi 


so  that  the  M  x  N  matrix  u>  of  the  un-weighted  averaging  example  is  given  by: 


(4.9) 


To  see  how  the  PSF  not  varying  spectrally  affects  IV,  we  partition  it  into  P  x  P  blocks 
Wij  where  1  <  i  <  M  and  1  <  j  <  N.  The  second  assumption  is  equivalent  to: 

Wlj=uijIP  (1  <i<M,  1  <  j  <  N), 

where  is  the  (i,j)  entry  of  uj  and  Ip  is  the  P  x  P  identity  matrix.  This  simple  structure 
of  Wij  will  be  important  when  the  Form  1  MAP  estimate  is  considered  in  conjunction  with 
principle  component  analysis. 

If  we  make  no  additional  assumptions  regarding  IF,  the  matrix  inversions  appearing 
in  the  MAP  estimates  are  potentially  full  MP  x  MP  matrices,  even  when  the  associated 
covariance  matrices  have  simple  structures  such  as  being  block  diagonal  with  P  x  P  blocks. 
In  order  to  avoid  working  with  full  MP  x  MP  matrices,  we  are  led  to  consider  a  third 
assumption: 

3.  The  PSF  combines  only  the  high  resolution  hyper- pixels  making  up  a  super-pixel  in 
order  to  form  the  corresponding  low  resolution  hyper-pixel. 
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This  is  equivalent  to  assuming  that  W  is  an  appropriately  generalized  version  of  the  form 
(4.6),  namely: 

M 

W  =  Q)Wm,  (4.10) 

m=  1 

where  Wm  is  the  P  x  LP  matrix  given  by: 

W  m  (^rra,(m— 1)L+1  Ipi  P  i  •  •  •  i  1  fC  TTI  fC  M.  (4.11) 

Here,  the  matrices  Wm  play  the  role  of  the  matrix  Wl  introduced  previously.  The  third 
assumption  is  characterized  by  the  condition  that  u  is  block  diagonal  with  1  x  L  blocks: 

^>m,k  7 ^  0  only  for  (m  —  1)L  +  1  <  k  <  mL  (1  <  m  <  M ), 

or  equivalently: 


u 


M 


where 


^ m  l)L+li  l)L+2j  ■  ■  •  i  ^m,mL)  • 

To  simplify  the  notation,  define  umj  =  com^m-i)L+j  so  that  u)  is  M  x  L  and 

^ m  j  &m,2i  ■  ■  •  1^1 n,L)  ? 

making  u T  the  mV1  row  of  d>. 

For  the  un- weighted  averaging  PSF  given  by  Equation  (4.9),  =  j-  for  1  <  rn  <  M 

and  1  <  j  <  L.  For  a  weighted  averaging  PSF,  we  may  define  the  M  x  L  matrix  Cj  in  any 
fashion. 

An  alternate  notation  for  Wm  is: 

Wm  = 

so  that 

M 

W  =  ©(i£®/p).  (4.12) 

m=  1 

Conditions  1,  2  &  3  are  summarized  by  the  assertion  that  each  entry  of  y  corresponds  to 
a  unique  band/super-pixel  pair  and  is  determined  from  y  =  Wz  +  n  as  a  weighted  average 
of  only  those  entries  of  z  that  lie  in  the  corresponding  band  and  super-pixel. 
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Throughout  the  remainder  of  this  chapter,  we  assume  all  three  of  the  above  conditions 
on  the  PSF.  If  they  don’t  hold  in  practice,  it  will  be  desirable  to  approximate  the  true  PSF 
with  one  where  they  do  hold. 

In  our  MATLAB  implementation  of  the  MAP  algorithm,  we  have  chosen  to  think  of 
a  PSF  as  a  Finite  Impulse  Response  (FIR)  filter.  Such  a  filter  may  be  represented  as  an 
Lv  x  Lh  matrix  where  Lv  =  nv/mv,  and  L/t  =  nh/mh  (so  L  =  LvLh).  In  the  present  context, 
this  is  equivalent  to  taking  all  rows  of  the  M  x  L  matrix  Co  to  be  identical,  where  each  row 
consisting  of  the  entries  of  the  Lv  x  Lh  FIR  matrix  representation  listed  in  column  order. 


4.4  Three  Forms  of  MAP  Estimators 


The  maximum  a  posteriori  (MAP)  estimate  (  [17],  p.350)  of  z  given  y  and  x  is  defined  by: 

z  =  argmaxPr  (z\xji) ,  (4-13) 

z 


where 


and  Pr  denotes  the  probability  density  function  of  the  continuous  random  variable  z\ip.  We 
derive  three  alternate  forms  of  this  estimate,  each  corresponding  to  it  own  function  of  z  to 
be  maximized.  It  will  be  seen  in  Section  4.5  that  all  three  are  equal  to  E(z\if>). 

The  three  forms  are  most  easily  derived  in  reverse  order  as  follows.  From  Bayes  rule: 


Pr{z\lp) 


PA^\z)Pr{z) 

PrW 


(4.14) 


where  we  have  chosen  not  give  different  names  to  different  density  functions  since  the  func¬ 
tion  argument  will  make  it  clear  which  PDF  is  being  referred  to.  Since  Pr  {xjo)  is  not  a 
function  of  2,  the  MAP  estimate  (4.13)  is  unaffected  if  it  is  ignored,  thus  giving  the  Form 
3  MAP  estimate: 


Zs  =  argmaxPr  (ij>\z)  Pr(z). 


(4.15) 


Making  the  reasonable  assumption  that  y\z  and  x\z  are  independent,  Equation  (4.14) 
yields: 


Pr(z  |  VO  = 


Pr  (y\z)  Pr  (x\z)  Pr{z) 


PrW 

Dropping  Pr  (ip)  as  before  yields  the  Form  2  MAP  estimate: 


(4.16) 


Z2  =  argmaxPr  (y\z)  Pr  (x\z)  Pr(z). 


(4.17) 
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Another  application  of  Bayes  rule  to  Equation  (4.16)  gives  the  following  representation: 


Pr(z\^)  = 


Pr  (y\z)  Pr  (z\x)  Pr(x) 


Pr  (</>) 

from  which  Pr  (i/>)  and  Pr(x )  may  be  dropped  to  obtain  the  Form  1  MAP  estimate: 


Z\  =  argmaxPr  (y\z)  Pr  (z\x) . 


(4.18) 


4.5  Derivation  of  Matrix  Equations 

From  Section  4.4,  the  Form  1  MAP  estimate  of  z  in  terms  of  probability  density  functions 
is  given  by: 


z  =  arg  max  Pr(y\z)Pr(z\x),  (4.19) 

z 

A  Bayesian  general  linear  model  is  assumed: 


y  =  Wz  +  n, 

meaning  that  the  noise  vector  n  is  independent  of  z  and  has  a  Gaussian  distribution  with 
mean  0  and  covariance  Cn.  This  leads  to  the  following  probability  density  function  for  y\z: 


Pr(y\z)  = 


:  exp 


--(y-wzyc-Hv-wz) 


(4.20) 


We  also  assume  the  Bayesian  general  linear  model  x  =  Sz  +  r/  which  implies  that  the 
joint  distribution  of  (x,  z)  is  multivariate  Gaussian.  By  Theorem  10.2  of  [17],  this  implies 
that  the  conditional  distribution  function  of  z\x  is  multivariate  Gaussian  as  well.  The 
relationship  between  the  joint  and  conditional  distributions  functions  is  presented  below. 

Let  yx  and  yz  denote  the  means  of  x  and  z  and  let  C(x.z)  denote  the  covariance  matrix 
of  ( x,z ).  Then, 


Pr(x,z ) 


1 

y/(2*)"+™\CM 


exp 


1 

2 


X  AG  A  q-i  f  x  P'x 
z-yz  )  PP  \  z-  nz 


5 


where  |G(X^)|  denotes  the  determinant  of  C(.JZZ) . 

The  covariance  matrix  C tXjZ\  is  an  (N  +  N P)  x  (N  +  NP)  real  symmetric  positive  dehnite 
matrix.  We  partition  C(X)Z)  as  follows: 


C 


(x,z) 


Cx 

Czx 


° zx 

Cz 


5 
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where  Cx  is  N  x  N,  Czx  is  NP  x  N,  and  Cz  is  NP  x  NP.  The  sub  matrices  Cx  and  Cz  are 
positive  definite  symmetric  since  C(X;Z)  is. 

By  Theorem  10.2  of  [17],  the  conditional  mean  vector  and  covariance  matrix  are  given 


by: 

/'.|,  =  +  CZzC~l(x  —  flj)  (4.21) 

C, „  =  C,  -  C^CpC^.  (4.22) 

From  the  model  x  =  Sz  +  r/,  we  have  the  following  relationships  (see  [17],  p.326): 

Hx  =  Syz  (4.23) 

a  =  SCzST4rCv  (4.24) 

czx  =  CZST,  (4.25) 

implying  that  the  conditional  mean  and  covariance  are: 

A h\x  =  Vz  +  CZST  [SCZST  +  Cv] _1  ( X  -  Sy,z)  (4.26) 

cz[x  =  Cz  -  CZST  [SCZST  +  cj_1  SCZ.  (4.27) 

Therefore,  the  conditional  PDF  of  z\x  is  given  by: 

1  1  T  1 

pr (z\x)  =  ^====exp^--(2:-^p)  C^(z-tiz  p)j,  (4.28) 

where  y,z\x  and  Cz\x  are  given  by  Equations  (4.26)  and  (4.27). 


It  is  seen  from  Equations  (4.20),  and  (4.28)  that  maximizing  the  product  Pr(y\z)Pr(z\x ) 
of  Equation  (4.19)  is  equivalent  to  minimizing  the  cost  function: 

Fi(z)  =  l(y-Wz)TC~\y-Wz)  +  ^(z-  yAx)T C~£(z-  nAx) ,  (4.29) 

so  that 

z  =  argminFKz).  (4.30) 

z 


The  cost  function  F\  may  be  minimized  by  the  standard  multivariate  calculus  technique 
of  taking  setting  the  gradient  of  F\  equal  to  zero  and  solving  for  critical  points.  Two  useful 
results  for  this  purpose  are  as  follows: 

g(z)  =  Uz~y)TMz~y)  =>  Vzd(z)  =  A(z~n) 

g(z)  =  \{y  —  W  z)T  A(y  —  W  z)  =>  Vzg(z)  =  -WTA(y-Wz). 


Applying  these  yield  the  following  Form  1  matrix  equation  whose  solution  is  the  value  of  z 
defined  in  Equation  (4.30): 


z  = 


wTc~1w  +  c;]l 


i-l  r 


wTc~ly  +  c;^zl 


-1 , 


(4.31) 
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The  matrix  inversion  lemma  (see  [17],  appendix  Al.1.3)  is: 


(A  +  BCD)-1  =  A"1  -  A~1B(DA~lB  +  C~1)-1DA~\  (4.32) 


where  the  indicated  inverses  are  assumed  to  exist  and  the  matrix  sizes  are  compatible  with 
the  indicated  products.  Applying  the  matrix  inversion  lemma  yields: 

£  =  [cz\x  -  CzlxWT(WCzlxWT  +  cn)-'wcz\x]  (wTc-1y  +  c;l1x^zlx). 

Let  E  =  WCZ\XWT  +  Cn  and  simplify  the  above  equation  as  follows: 


Ai 

X 

Cz 

xwte~ 

Hvc 

z\x] 

(wTc-'ypc;]l 

Pz\x) 

Vz\x 

+ 

Cz\x 

WT 

[c- 

ly- 

E~ 

1WCz\xWTC~1y 

-E~ 

lWyz 

\x ] 

l1 z\x 

+ 

Cz\x 

WT 

[c„- 

ly  - 

E~ 

'(B-CyC-1  y 

-  E~ 

1Wy,z 

J 

Fz\x 

+ 

Cz\x 

WT 

[c„- 

ly  - 

(' 

-  E-'Cn)  C~ly  - 

-E~a 

■w^ 

J 

A6 z\x 

+ 

Cz\x 

WT 

[Q 

ly  - 

c- 

‘j/  +  E-'y  -  E - 

lWyz 

'\x\ 

ft'zlx 

+ 

Cz\x 

WT 

E~' 

( y - 

Therefore, 

*  =  +  C,„Wt[WCm,Wt  +  Cn}-\ y  -  (4.33) 

which  is  the  starting  point  of  Section  4.11. 

From  Section  4.4,  the  Form  2  MAP  estimator  in  terms  of  PDF  functions  is  given  by: 


z  =  &vgm&xPr(yy\z)Pr(x\z)Pr(yz), 
z 


where 


Pr{x\z)  = 
PT(z)  = 


i 


exp 


1 

2* 


exp 


)Tc-\ 

\T  r*— 1 


Z~Pz)  Cz  (Z-Pz) 


(4.34) 

Sz)  , 

(4.35) 

Pz)  ■ 

(4.36) 

Analogous  to  Form  1,  the  cost  function  for  Form  2  is  obtained  from  Equations  (4.34),  (4.20), 
(4.35),  and  (4.36): 


=  \(v-Wz)TC-\y-Wz)+l-(x-Sz)TC-\x-Sz) 


1 


+  ^(z-  yz)T  Czx{z-  yz). 


(4.37) 


Minimizing  the  cost  function  F2  of  Equation  (4.37)  yields  the  Form  2  MAP  estimate: 

z  =  [c;1  +  wTc-1w  +  sTc-1s]~1[c;1nz  +  wTc-1y  +  sTc~1x],  (4.38) 
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which  is  identical  to  Equation  (A.l)  and  is  the  starting  point  of  Appendix  A.  In  Appendix 
A,  the  matrix  inversion  lemma  is  applied  twice  in  succession  resulting  in  a  simplified  version 
of  the  above  Form  2  matrix  equation. 

From  Section  4.4,  the  Form  3  MAP  estimator  in  terms  of  PDF  functions  is  given  by: 


z 


argmaxPr  (ip\ z)  Pr(z), 
z 


where  ip 


Making  the  associations  y  * — >  ip  and  z\x  < — >  z,  we  see  from  Equation  (4.19)  that  Form 
3  and  Form  1  have  the  same  structure.  If  we  further  associate: 


W 

cn 


( I ) 

(  cn  0  \ 

1  0  cj’ 


then  we  may  use  Equation  (4.29)  to  obtain  the  Form  3  cost  function: 


* «*>  =  \ 


y 

X 


W 

s 


Cn  0 

0  Cn 


-1 


y 

X 


W 

s 


+  ^(z-  Vzf  Czl(z-  nz), 


and  use  Equation  (4.33)  to  obtain  the  Form  3  matrix  equation: 


z  —  +  Cz 


We  have  assumed  that  17  and  n  are  independent  in  order  to  use 
ance  matrix  of  ip\z. 


Cn 

0 


(4.39) 


as  the  covari- 


Since  Form  1  is  the  only  form  that  is  applicable  when  the  model  x  =  S z  +  77  is  not 
assumed,  the  following  theorem  proves  that  it  is  sufficient  to  program  just  Form  1  in  order 
to  incorporate  all  three  forms  of  the  MAP  algorithm. 


Theorem  4.5.1  Under  the  assumption  x  =  Sz  +  77,  the  three  forms  of  MAP  estimators 
defined  above  are  equivalent  and  each  is  equal  to  E(z\ip). 

Proof  The  Form  2  and  Form  3  MAP  estimates  are  easily  seen  to  be  identical  by  checking 
that  their  associated  cost  functions  F^iiz)  and  Ffiz)  are  identical.  We  may  prove  that  the 
Form  1  and  Form  2  MAP  estimates  are  equivalent  either  by  proving  that  the  gradient  of  the 
difference  of  their  cost  functions  is  zero,  or  by  proving  the  equivalence  of  Equations  (4.31) 
and  (4.38)  directly.  The  latter  approach  is  taken  below. 
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Applying  the  matrix  inversion  lemma  to  Equation  (4.27)  yields: 

1  _  1  ,  nT/i-1  n 

° z\x  —  -h  0  0^0, 

and  substituting  this  into  Equation  (4.31)  gives: 

z  =  [W'c-^w  +  C71  +  STC-1S]  ■'  [iWCpy  +  (Cp  +  SU^'S)  . 

Comparing  the  above  with  Equation  (4.38),  we  see  that  the  Form  1  and  Form  2  MAP 
estimates  are  equivalent  provided  that: 

(cy  +  src^s)  m*  =  +  sTc,-u, 

where  /nz|x  is  given  by  Equation  (4.26).  The  following  calculation  verifies  this: 

(C71  +  sTc~1s)  pzV x  =  (cy  +  sTc~ls )  [»z  +  czsT  [sczsT  +  cv] _1  (*  -  s»z) 

=  C~1iiz  +  ST  [ SCZST  +  Cv] _1  (*  -  Aju  J  +  STC~lS^z 
+  SFC^SCtSf1,  [SCZST  +  CVYX  (x  -  A/xJ 
=  C~l[iz  +  AT  [AC2At  +  CJ _1  (*  -  A/x;)  +  A'C,,  '5//, 

+  STC~l  [(SCZST  +  Cv)  -  Cv]  [SCZST  +  Cv] _1  (x  -  Snz ) 

=  C- V,  +  AT  [AC;At  +  C„] _1  (®  -  A/*;)  +  STC~lSnz 
+  [s^C"1  -  [ACzAr  +  Cj  _1]  (aj  -  A/xJ 

=  c;yz  +  sTcyx. 

Since  we  have  shown  that  all  three  forms  are  identical,  the  fact  that  each  is  equal  to 
E(z \ij))  is  established  by  showing  it  for  any  one  of  the  forms.  A  straightforward  application 
of  Theorem  10.3  of  Kay  [17]  shows  that  the  Form  3  MAP  estimate  given  by  Equation  (4.39) 
is  equal  to  E(z \ip).  m 


4.6  Estimation  of  the  Spectral  Response  Matrix 


In  Section  4.7,  we  will  use  the  relationship  x  =  Sz  +  r/  to  estimate  the  covariance  matrix 
Czix  in  terms  of  the  covariance  matrix  Cz.  For  this,  it  will  be  important  to  know  how  the 
spectral  response  matrix  A  is  estimated  and  what  its  resulting  structure  is. 


For  each  low  resolution  hyper-pixel  yrn ,  let  xm  denote  the  average  of  the  corresponding 
broadband  image  data: 


Xm=L 


xq(( 


m—l)L+j)  • 


(4.40) 
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Consider  the  equations: 


Xm  =  ( Vrn , s)  (1  <  m  <  M),  (4.41) 

where  the  angle  brackets  denote  the  standard  inner  product  and  s  =  (si,S2,  •  •  •  ,  s.p)T.  In 
matrix  form  Equation  (4.41)  is: 

Ys  =  x,  (4.42) 

where  x  =  (xi,x2, ,  xm)t  and  Y  is  the  M  x  P  matrix  determined  by  placing  yfn  in  row 
m  for  1  <  m  <  M. 

Equation  (4.42)  is  solved  in  the  least  squares  sense,  subject  to  i  sj  =  I  and  sj  Y  0  for 
1  <  j  <  P,  in  order  to  determine  a  normalized  spectral  response  vector  s.  The  normalized 
spectral  response  vector  is  used  to  form  an  N  x  NP  matrix  S  =  SQT  as  follows: 

N 

S  =  0sT  =  I/v  ®  sT.  (4.43) 

n= 1 

Observe  that  for  the  no  noise  case,  x  =  Sz  +  77  is  equivalent  to: 

Xi  =  (zi,  s)  (1  <  i  <  N ), 

which  may  be  viewed  as  Equation  (4.41)  extended  from  low  resolution  hyper-pixels  to  high 
resolution  hyper-pixels.  Also  observe  that  the  assumed  form  of  S  dictates  that  aq  is  de¬ 
termined  solely  from  r/,;  and  the  hyper  pixel  zL  corresponding  to  xt.  Thirdly,  note  that 
x  =  Sz  +  r)  determines  x  t]  in  super-pixel  ordering  while  x  =  Sz  +  77  determines  it  in 
lexicographical  order. 

In  MATLAB,  the  computation  of  s  is  given  by: 

s  =  LSQLIN(Y,  x,  [  ],[  ] ,  ones(l,  P),  1,  zeros(l,  P),  ones(l,  P)); 

where  LSQLIN  is  a  function  available  in  the  MATLAB  optimization  toolbox. 

SlHC6 

||Ts  —  x\\2  =  stYtYs  —  2xtYs  +  xTx, 

the  problem  may  be  formulated  as  a  convex  quadratic  optimization  problem.  In  this  for¬ 
mulation,  let  C  =  YTY,  d  =  —YTx  and  compute  s  such  that: 

^sr(7s  +  dTs 

_ ^  P 

is  minimized  subject  to  the  linear  constraint  l2j=1  sj  =  1  and  the  bound  constraints  Sj  >  0, 
1  <  j  <  P.  The  freely  available  C  program  QLD  [20]  may  be  used  to  solve  this  equivalent 
problem.  We  have  determined  through  numerical  experimentation  that  the  computations 
of  C  and  d  incur  excessive  round  off  errors  if  performed  in  single  precision. 
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4.7  Estimation  of  Conditional  Covariance  from  Un¬ 
conditional  Covariance 

In  this  section,  we  assume  the  linear  model  x  =  Sx  +  rj  holds,  which  implies  that  the 
relationship  between  the  covariance  matrices  Cz\x  and  Cz  given  by  Equation  (4.27)  holds: 

cz\x  =  cz-  CZST(SCZST  +  CJ-'SC,.  (4.44) 

Therefore,  the  relationship  between  the  covariance  matrices  Cz\x  and  C~  is  given  by: 

Cz,x  =  C~z  -  CZST(SCZST  +  C^SCs  ,  (4.45) 

where  S  =  SQT  ■  Cz  =  QCZQT ,  Cz\x  =  QCZ\XQT ,  and  z  gives  the  super-pixel  ordering  of  z 
as  defined  in  Section  4.2. 

An  assumption  that  will  facilitate  the  parallel  implementation  of  the  Form  1  MAP 
estimate  is  that  Cz \x  be  block  diagonal  with  P  x  P  blocks.  Due  to  the  special  form  of  S 
given  by  Equation  (4.43),  we  will  see  below  that  Cz\x  has  this  form  if  C~  has  it.  Accordingly, 


N 

C4  =  0fit,  (4.46) 

2=1 

where  each  Bi  is  a  P  x  P  matrix.  From  Equation  (4.43),  the  N  x  N  matrix  SCzST  is  given 
by: 

N 

SCZST  =  0stB,s. 

2=1 

Assuming  that  Cv  =  it  follows  that 

N 

(SCZST  +  C,)-1  =  Q)d-\ 

2=1 

where  dt  =  (sTBjS  +  a^).  Therefore  Equation  (4.44)  becomes 

c-z\x  =  cz-  CZSTP>~XSCZ, 

where  D  =  d*.  Since  STD~1S  =  d^s sT,  it  follows  from  Equation  (4.46)  that: 

N 

cz\x  =  0  [Bi  -  dr'BiS  (5iS)T]  ,  (4.47) 

2=1 

and  we  see  that  the  desired  structure  for  Cz\x  is  inherited  from  Cz . 
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Starting  with  Equation  (4.26),  the  following  expression  for  is  obtained  in  a  similar 
way: 

_  (  T  T  T  \T 

\x  \l^Zl\x^Z2\x^  '  '  '  ^ZN\x)  -) 

Hzi\x  =  Hi  +  (xf  -  s T/ui)d~1Bis,  (4.48) 

where  p*  is  a  P  x  1  vector  that  denotes  the  mean  of  the  high  resolution  hyper-pixel  A  =  zq^ . 
Note  that  the  same  coefficient  d^BiS  appears  in  both  (4.47)  and  (4.48)  which  may  be 
utilized  for  computational  efficiency. 

For  1  <  m  <  M  and  1  <  j  <  L,  define  the  P  x  P  matrix  Bm  j  =  The  Bmj 

merely  represent  an  alternative  indexing  of  the  Bi,  convenient  when  identifying  super-pixel 
components.  Thus, 

M 

Cz  =  Bm,  (4.49) 

m= 1 
L 

Bm  =  (4.50) 

3= 1 

A  similar  notation  is  used  for  <77^1;,;: 

M 

Cz\x  =  0Gm)  (4-51) 

m= 1 
L 

Gm  =  (4.52) 

3= 1 

Equation  (4.47)  implies: 

Gm,j  —  Bmj  —  (sT BmjS  +  cr^)  1  (Bmj s)  (Bmj s)T  .  (4.53) 

Starting  with  estimates  Bmj  we  use  Equation  (4.53)  to  obtain  Gmj  which  results  in  an 
estimation  of  Cz\x-  Thus,  obtaining  a  good  estimate  of  C%\x  comes  down  to  obtaining  good 
estimates  of  the  Bmj.  There  is  a  caveat  however:  the  matrix  given  by  Equation  (4.53)  is 
singular  when  av  =  0.  To  see  this,  let  am  ]  =  sTBmjS  and  v  =  Bmjs.  It  is  easily  verified  that 
amj  is  an  eigenvector  of  Bmj ssT  with  eigenvector  v.  Therefore,  det (amjlp  —  £>mjssT)  =  0, 
which  implies  that  Gmj  is  singular  when  an  =  0. 

The  Form  1,  2,  and  3  MAP  estimates  are  sensitive  to  a  condition  we  refer  to  as  the 
special  case  condition  (or  assumption).  Explicitly,  this  is  the  condition: 

Bm,j  =  Bm, i  (1  <  m  <  M,  1  <  j  <  L). 

We  colloquially  refer  to  this  by  stating  that  the  P  x  P  diagonal  blocks  of  C7  are  identical 
within  super-pixels.  It  is  readily  seen  from  Equation  (4.53)  that  Cz\x  inherits  the  special 
case  condition  from  67 . 
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4.8  Consequences  of  Making  the  Conditional  Covari¬ 
ance  Matrix  Diagonal 

In  this  section  we  discuss  how  making  Cz\x  a  diagonal  matrix  implies  an  undesired  restric¬ 
tion  on  the  spectral  response  matrix  S.  In  particular,  such  a  diagonal  structure  imposes 
restrictions  that  are  most  likely  inconsistent  with  an  estimation  of  S  such  as  that  given  in 
Section  4.6. 

We  begin  by  examining  Equation  (4.53)  for  the  case  that  Brn  j  is  a  diagonal  matrix. 
In  this  case,  it  is  shown  below  (see  Theorem  4.8.1)  that  s  has  only  one  nonzero  entry, 
which  must  be  1  due  to  the  normalization  Ylj=i  sj  =  1-  As  another  case  (see  Theorem 
4.8.2),  assuming  (4.53)  is  diagonal  when  P  =  2  leads  to  the  result  that  SiS2  is  a  multiple 
of  cA  where  the  multiple  depends  only  on  the  entries  of  Bmj.  Either  case  leads  to  the 
conclusion  that  Gni  J  (and  hence  the  conditional  covariance  matrix  C~\x)  may  not  be  taken 
as  a  diagonal  matrix  without  limiting  the  entries  of  s  in  a  manner  most  likely  inconsistent 
with  an  estimation  of  S  such  as  that  given  in  Section  4.6. 


Theorem  4.8.1  IfGm  j  and  Brn  j  are  diagonal  then  s  =  e*  for  some  i  e  {1,  2, . . . ,  P}  where 
e*  denotes  the  unit  vector  with  a  one  in  its  i^1  entry  and  zeros  elsewhere. 

Proof  From  Equation  (4.53),  we  see  that  Bmjs(Bmjs)T  is  diagonal.  Let  v  =  Bmj s  so 
that  vnVk  =  0  for  n  ^  k.  Since  s  ^  0,  v  ^  0  and  we  may  choose  i  such  that  u,  ^  0.  There¬ 
fore,  v  =  rye,.  However,  v  =  (biSi,b2s2, . . . ,  bpSp)T  where  £>mj  =  diag(&i,  b2, . . . ,  bp),  so 
it  follows  that  bkSk  =  0  whenever  k=fi.  Since  Bmj  is  positive  definite,  /y  >  0  for  all  k. 
Therefore,  Sk  =  0  for  k  ^  i  and  since  JT=1  Sj  =  1,  s;  =  1.  ■ 


Theorem  4.8.2  If  Gmj  is  diagonal  and  P  =  2,  then 

SlS2  =  ~W 

where 


Proof  Using  Equation  (4.53),  we  set  the  (1,  2)  entry  of  Gmj  equal  to  zero: 

^1,2  (^1,1^1  T  2&i)2SiS2  T  b2  2S2  T  c^^)  1(^i,i,5i  T  bi,2s2)(bi,2Si  T  b2,2s2)  =  0. 

and  the  result  easily  follows.  ■ 

For  the  general  case  when  Bm  ]  is  any  positive  definite  symmetric  matrix  and  P  is  any 
positive  integer,  we  argue  heuristically  as  follows.  Setting  the  off  diagonal  entries  of  the 
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symmetric  matrix  Gmj  equal  to  zero  leads  to  the  following  system  of  P(P  —  l)/2  nonlinear 
equations  in  P  unknowns: 

p  p 

^  ^  ~^(bj,nbk,£  bpnb^j£)SgSn  ^p^i,ki  \  ^  i  k  P. 

£=1  n= 1 

Because  of  the  condition  Xq=i  sj  =  1>  each  °f  these  equations  represents  a.  P—2  dimensional 
surface  that  lies  in  a  P  —  1  dimensional  space.  The  intersection  of  all  P(P  —  l)/2  surfaces 
represents  the  set  of  solutions  s  for  which  Gmj  is  diagonal.  Intersecting  these  surfaces  will 
result  in  a  solution  set  of  dimension  P  —  2  or  less,  most  likely  of  dimension  much  less  than 
P  —  2  for  sufficiently  large  P.  (We  see  from  the  examples  given  above  that  solution  spaces 
of  dimension  0  occur  when  Bmj  is  diagonal  or  when  P  —  2.)  We  conclude  that  since  the 
solution  set  is  of  dimension  P  —  2  or  less,  a  procedure  for  determining  s  such  as  that  of 
Section  4.6  will  most  likely  be  inconsistent  with  the  assumption  that  conditional  covariance 
matrix  C%\x  is  diagonal. 


4.9  Conditional  Independence 

This  section  defines  two  levels  of  conditional  independence  and  relates  them  to  the  assumed 
structures  of  associated  matrices.  The  first  level  factors  the  probability  distribution  function 
of  z\x  into  the  conditional  distribution  functions  of  zn \ xn ,  where  zn  is  the  nl  hyper-pixel 
of  z,  and  the  second  level  factors  it  further  into  the  conditional  distribution  functions  of 
zp,n\xn  where  zn  =  (zi^n,  z^n,  •  •  • ,  zP,n)T ■  The  first  level  of  conditional  independence  follows 
from  assumptions  previously  made  concerning  the  P  x  P  diagonal  block  structure  of  the 
covariance  matrix  Cz  and  the  Bayesian  general  linear  model  x  =  Sz  +  rj,  while  the  second 
level  is  equivalent  to  the  assumption  that  the  conditional  covariance  matrix  Cz\x  is  diagonal. 


Explicitly,  the  two  levels  of  conditional  independence  are: 

Level  1  :  Pr(z\x ) 

N 

J~J  ^>r  fan \%n)  i 

71=1 

(4.54) 

Level  2  :  Pr(zn \xn) 

II 

5 

s 

H 

3 

II 

to 

.,N. 

(4.55) 

We  show  below  that  conditional  independence  level  1  is  implied  by  the  Bayesian  general 
linear  model  x  =  Sz  +  t]  and  the  assumed  structures  of  Ch  and  S.  Since  Cz  is  block  diagonal 
with  P  x  P  blocks,  and  the  entries  of  Cz  corresponding  to  zn  and  z^  are  proportional  to 
the  correlations  p(zitn,  Zjtk)  where  1  <  i.j  <  P  and  1  <  n,  k  <  N,  it  follows  that  these 
correlations  are  zero  whenever  n  ^  k.  We  express  this  succinctly  as:  p(zn,  =  0  whenever 
(n  7^  k). 
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From  Equation  (4.24)  and  Cz  =  QrC~Q,  it  easily  follows  that: 

cx  =  SCZST  +  Cv, 

where  S  =  SQT  is  given  by  Equation  (4.43).  This  and  Equation  (4.46)  imply  that  Cx  is  a 
diagonal  matrix: 

N 

cx  =  ®  (s' TBns)  +  C,r 

n= 1 

so  that  p(xn,xk )  =  0  whenever  n  ^  k.  Similarly,  equation  (4.25)  easily  implies: 

/  N 

Czx  =  QT  I  Bns 

\n= 1 

so  that  p(xn,  zk)  =  0  whenever  n  /  k.  Collecting  the  above  results  we  have: 

p(xn,xk)  =  p(xn,zk)  =  p(zn,zk )  =  0  (n^k). 


For  Gaussian  distributions,  a  zero  correlation  between  two  random  variables  implies 


independence  between  them.  Therefore, 

Pr(x)  = 

N 

Y[Pr(xn) 

- 1 

II 

TT 

it - 1 

N 

Pr  (xn,  £n) 

n= 1 

and  equation  (4.54)  readily  follows: 

Pr(z\x)  = 

Pr{x,z) 

Pr{x) 

TT  Pr(xn,  Zn) 

n=l  Pr(Xn) 

N 

J  J  Pr  \%n)  • 

n= 1 

Since  z\  and  x\  are  independent  of  (x2,  x3, . . . ,  xjy)  we  have: 


Pr(zi\x) 


PyjZuX) 

Pr(x )  ’ 

Pr(z1,X1)Pr(x2,X3,  ...,XN) 
Pr(x1)Pr(x2,x3, .  .  .  ,XN) 

Pr(zi\xi), 
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and  similarly, 


Pr{,^n\'^)  Pr(zn\xn),  (1  ^  77.  ^  N^j . 

Arguing  in  the  same  manner  yields: 

Pr (Zp,n \x)  Pr (zp,n \xn) >  (1  E  72  E  A^,  1  5;  P  S  T*)- 

Using  the  fact  that  for  Gaussian  distributions,  zero  correlation  and  independence  are 
equivalent,  it  follows  that  Cqx  is  diagonal  if  and  only  if 

p 

Pr  (zn  |  *c)  1 1  Pr  i.zp,n  |*®)  >  ^  1)2,...,  N. 

P= 1 

From  the  above  results  it  follows  that  C~\x  is  diagonal  if  and  only  if  conditional  independence 
level  2  (i.e.  Equation  (4.55))  holds  for  n  —  1, 2, . . . ,  N . 

As  remarked  in  Section  4.8,  the  assumption  that  Cx\x  is  diagonal  is  likely  to  be  inconsis¬ 
tent  with  the  estimation  given  in  Section  4.6,  making  level  2  conditional  independence  an 
undesirable  assumption. 


4.10  Principal  Component  Analysis 

Principal  component  analysis  (PC A)  starts  with  the  construction  of  a  P  x  P  “spectral” 
covariance  matrix  C  that  represents  the  covariance  among  the  entries  of  a  high  resolution 
hyper-pixel.  Since  z  is  presumed  unknown,  the  hyper-pixels  of  y  may  be  used  to  compute 
an  estimate  of  C.  Using  an  eigenvalue  decomposition  algorithm,  an  orthonormal  set  of 
eigenvectors  of  C  is  obtained  and  is  used  to  construct  an  orthogonal  matrix  E  whose  columns 
are  the  eigenvectors.  The  unknown  high  resolution  hyper-pixels  z,  (i  =  1,2,  ...,N)  are 
related  to  a  different  set  of  unknowns  z,-  through  the  transformation  z,  =  Ezt.  The  problem 
of  computing  z  such  that  y  =  lUz  +  n  and  x  =  Sz  + 1]  is  transformed  into  computing  z 
such  that  y  =  Wz  +  n  and  x  =  Sz  +  rj,  where  the  quantities  with  bars  are  defined  below. 

Let  En  =  0^=1  E  and  EM  =  ®^=i  E  and  define  z  =  Ej^z ,  y  =  E^y,  n  =  E^ n, 
S  =  SEn,  and  tU  =  EJjWEn-  Since  E~x  =  ET ,  the  same  is  true  of  En  and  Em,  making 
it  straightforward  to  show  the  equivalence  of  x  =  Sz  +  t]  and  x  =  Sz  +  y  as  well  as  the 
equivalence  of  y  =  TUz  +  n  and  y  =  i-Uz  +  it.  Since  n  represents  white  Gaussian  noise, 
so  does  n  and  it  follows  that  an  image  enhancement  algorithm  for  computing  z  may  be 
applied  to  obtain  z,  where  the  relevant  input  quantities  are  substituted  with  their  barred 
counterparts.  We  also  observe  that  since  E  is  orthogonal,  ||h||  =  ||n||  so  that  noise  is 
not  magnified  when  working  in  the  principal  component  domain.  Although  z  may  be  used 
directly  for  producing  images,  z  may  be  obtained  from  z  =  ENz  if  desired. 
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The  spectral  response  matrices  in  their  respective  domains  being  related  by  S  =  SEn 
implies  that  the  P  x  1  corresponding  spectral  response  vectors  are  related  by: 

s  =  Ets  (4.56) 

If  the  PSF  matrix  W  satisfies  assumption  2  of  Section  4.3  (i.e.  does  not  vary  spectrally), 
then  W  =  W.  This  is  important  since  many  of  the  image  enhancement  algorithms  we 
have  examined  take  advantage  of  special  assumed  forms  of  W  for  efficiency  purposes.  In 
particular,  the  current  version  of  the  Form  1  MAP  algorithm  makes  use  of  PSF  assumptions 
1,  2,  and  3  of  Section  4.3  for  efficient  coding.  We  prove  below  that  assumption  2  implies 

W  =  W. 

As  in  Section  4.3,  we  partition  W  into  P  x  P  blocks  Wij  where  1  <  i  <  M  and 
1  <  j  <  N.  It  is  easily  seen  that  the  corresponding  (i ,  j )  block  of  W  is  given  by  ETWt  JE. 
The  assumption  that  W  does  not  vary  spectrally  is  equivalent  to  Wij  being  a  multiple  of  Ip 
for  all  Therefore,  WhJ  commutes  with  E.  Since  E  is  orthogonal,  ETE  =  Ip.  Hence, 

the  (i,j)  block  of  W  is  implying  that  W  =  W. 

Since  the  eigenvectors  of  E  corresponding  to  the  smaller  eigenvalues  are  not  expected  to 
contribute  to  the  enhancement  of  images,  we  now  examine  a  “limited”  PC  A  where  only  the 
largest  eigenvalues  are  computed.  The  intent  is  to  improve  performance  without  degrading 
image  quality.  Accordingly,  choose  p  (1  <  p  <  P)  and  compute  the  eigenvectors  of  C  belong¬ 
ing  to  the  largest  p  eigenvectors.  Let  E  be  the  P  xp  matrix  defined  by  E  =  (Ei,  E2, ,  Ep) 
where  the  Ei  are  P  x  1  eigenvectors  sorted  by  eigenvalue  from  largest  to  smallest.  Using 
the  same  definitions  given  in  above,  E ^  becomes  NP  x  Np,  Em  is  MP  x  Mp,  z  is  Np  x  1, 
S  is  N  x  Np,  y  is  Mp  x  1,  and  n  is  Mp  x  1. 

After  applying  an  estimate  such  as  the  Form  1  MAP  estimate,  we  obtain  an  estimate 
of  z.  As  before,  z  may  be  used  directly  to  obtain  images.  Optionally,  the  corresponding 
estimate  of  z  is  given  by  z  —  E^z.  We  note  that  the  resulting  size  of  z  is  NP  x  1  even  if 
p  <  P.  For  efficiency,  the  product  EN z  may  be  obtained  from  the  product  E[z\,  z%,  ■  ■  ■ ,  Av], 
where  each  z*  is  p  x  1  and  z  =  {zp  z2\  ■  ■  ■ ;  Zn)-  The  resulting  P  x  N  matrix  may  then  be 
reshaped  to  a  NP  x  1  vector  to  obtain  z. 

There  is  a  less  efficient,  but  perhaps  more  intuitive  way  to  perform  limited  PC  A  analysis. 
In  this  formulation,  y  is  kept  as  MP  x  1,  but  with  each  low  resolution  hyper-pixel  having 
a  zero  in  its  last  P  —  p  entries.  We  also  leave  E  as  a  P  x  P  matrix  but  take  its  last  P  —  p 
columns  to  be  zero  vectors  of  length  P.  This  results  in  essentially  the  same  z  as  above,  but 
now  each  Zj  is  a  P  x  1  vector  with  zeros  in  the  last  P  —  p  entries.  The  resulting  limited 
PCA  computation  of  z  is  identical  to  the  one  obtained  above.  The  advantages  of  the  former 
method  is  not  only  that  it  takes  less  storage,  but  that  it  also  takes  less  time  to  compute, 
because  calculations  that  multiply  zero  times  itself  are  eliminated. 

The  advantage  of  PCA  space  processing  lies  in  the  fact  that  most  of  the  signal  energy  lies 
in  the  PCA  subspace  corresponding  to  the  top  eigenvalues.  The  lower  principal  component 
dimensions  can  be  treated  with  a  simpler  algorithm,  such  as  interpolation,  or  even  set  to  zero 
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if  they  are  deemed  to  be  primarily  noise.  Processing  in  a  lower  dimension  PCA  subspace 
reduces  the  size  of  the  matrix  inverse  required  in  the  closed  form  solution  and  requires  the 
estimation  of  fewer  statistical  parameters.  In  addition  to  the  benefit  in  terms  of  processing 
time,  this  tends  to  lead  to  covariance  estimates  that  are  better  conditioned,  making  the 
resulting  MAP  algorithm  more  robust. 


4.11  Analysis  of  the  Form  1  MAP  Estimate 

Starting  with  the  Form  1  MAP  estimate  presented  in  Section  4.5,  this  section  derives  an 
equivalent  expression  that  is  suitable  for  parallel  implementation.  Assuming  the  special 
case  where  the  P  x  P  blocks  of  CP,  are  identical  within  super-pixels,  we  then  prove  that  the 
matrix  appearing  as  an  inverse  in  the  new  expression  is  singular  if  and  only  if  an  =  av  =  0 
(i.e.  there  is  no  noise).  A  separate  expression  is  developed  under  the  noise  free  and  special 
case  assumptions  and  in  this  expression,  no  matrix  inverses  whatsoever  appear.  Having 
no  matrix  inverses  to  compute  leads  to  an  algorithm  whose  complexity  is  low  enough  to  be 
considered  for  real  time  implementation  on  a  parallel  computer  system. 

In  order  to  facilitate  parallel  computations,  each  of  these  Form  1  MAP  expressions  are 
broken  down  into  super-pixel  components.  This  allows  the  work  to  be  distributed  among 
the  available  processors  in  a  convenient  and  efficient  fashion.  Throughout  this  section,  the 
three  conditions  of  Section  4.3  are  assumed  and  the  definitions  of  Bm  j  given  by  Equations 
(4.49)  and  (4.50)  are  adopted. 

As  derived  in  Section  4.5,  the  Form  1  MAP  estimator  is  given  by: 

z  =  nAx  +  CzlxWT[WCzlxWT  +  Cn\~\y  -  Wpizlx) 

If 

Z  =  Qz ,  Hz\x  =  QHz\xi  C%\x  =  QCZ\XQ  i 

where  Q  is  the  super-pixel  permutation  matrix  defined  in  Section  4.2,  then  Equation  (4.57) 
becomes 

z  =  n-Ax  +  Cz\xWT[WC-z\xWT  +  Cn]~\y  -  WixAx),  (4.58) 

and  we  see  that  Equations  (4.57)  and  (4.58)  have  the  same  form.  Equation  (4.58)  is  much 
easier  work  with  than  (4.57)  due  to  the  simple  structure  of  W.  The  estimate  z  is  given  by 
z  =  QTz  which  may  be  implemented  in  an  efficient  manner  by  permuting  the  hyper-pixels 
of  z  as  follows: 


(4.57) 


%  % q~1(2)t  •  ■  ■  i  Zq-l(N))  > 

where  the  semi-colon  denotes  vertical  concatenation  and  q  denotes  the  super-pixel  permuta¬ 
tion  corresponding  to  Q  defined  in  Section  4.2.  In  the  notation  above,  zg- qn)  =  zn  represents 
tb 

the  nl  hyper-pixel  for  n  =  1,  2, . . .  N . 
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A  simplifying  assumption,  which  appears  to  be  necessary  for  efficient  parailei  compu¬ 
tations,  is  that  Cz  =  QCZQT  be  block  diagonal  with  P  x  P  blocks.  From  Section  4.7,  we 
know  that  Cz  having  this  structure  implies  that  Cz\x  =  QCZ\XQT  will  as  well.  Since  pre 
multiplication  by  Q  followed  by  post  multiplication  by  QT  permutes  the  P  x  P  blocks  along 
the  macro  diagonal  of  such  matrices,  it  follows  that  Cz ,  Cz\x ,  and  Cz\x  all  inherit  the  same 
P  x  P  block  diagonal  structure  from  Cz. 

Let  L  =  N/M  and  partition  Cz\x  into  diagonal  blocks  Gi,  G2,  ■ . . ,  Gm  where  each  Gm  is 
a  LP  x  LP  block  diagonal  matrix  with  P  x  P  blocks.  Thus, 

M 

Cz\x  =  0Gm,  (4.59) 

m= 1 

and  each  Gm  for  m  =  1,  2, . . . ,  M  has  the  form: 

L 

Gm  =  (^Grn:j,  (4.60) 

5= 1 

where  Gmj  is  a  P  x  P  matrix.  It  follows  from  Equations  (4.12)  and  (4.59)  that 

M 

WCzlxWT  =  0  Gm 

m= 1 

where  Gm  is  the  P  x  P  matrix  given  by: 

Gm  =  ®  Ip)Gm(um  ®  Ip)  =  (4-61) 

Due  to  the  simplifying  assumption  that  the  noise  covariance  matrix  Cn  is  a  multiple 
of  the  MP  x  MP  identity  matrix,  the  inverted  matrix  of  Equation  (4.58)  becomes: 

M 

[WC-Z\XWT  +  Cn]-1  =  ®  (Gm  +  a2nIp)-\  (4.62) 

m=  1 

Similarly,  the  NP  x  MP  product  CZ\XWT  appearing  in  Equation  (4.58)  is  given  by: 

M 

CZ\XWT  =  Gm(ujm  ®  Ip).  (4.63) 

m— 1 

The  matrix  CZ\XWT  of  Equation  (4.63)  is  a  non  square  (NP  x  MP)  matrix  that  consists  of 
M  blocks  on  the  diagonal,  each  of  size  LP  x  P. 

Using  the  superscript  (1)  to  designate  a  quantity  that  pertains  to  the  Form  1  MAP 
estimate,  define  the  NP  x  MP  matrix  by: 


so  that  the  Form  1  MAP  estimate  is  given  by: 


z  =  nz 


+  Am(y-W^,x). 


Putting  Equations  (4.62)  and  (4.63)  together  yield: 


Gm(u}m  0  Ip)(Gm  +  cr^Ip) 


®A (!) 


The  LP  x  P  matrix  Am  is  defined  by: 


^■rn  —  Ip)(Gm  +  cr^lp) 


which  when  written  out  becomes: 


aL1}  = 


&m,lGm,l(Gm  +  anlp)  1 
+  an^p)  1 

& m,L G  m,L  (G  m  +  (T^Ip)  1 


From  Equation  (4.58),  we  see  that  the  Form  1  estimator  is  given  by: 


®  4J1  (y  -  |, 


2  =  AA 


Breaking  this  up  into  M  super-pixel  components  yields: 


or  equivalently: 


Pzm\x  +  Am  [Vm  -  {Am  ®  Ip) 


=  [hp  ~  Am  (^m  ®  Ip)]  Vzm\x  +  AmVm, 


(4.64) 


(4.65) 


(4.66) 


(4.67) 


where  zm  and  fj,Bm \x  are  LP  x  1  and  ym  is  a  P  x  1  hyper-pixel.  The  breakdown  of  z  into 
M  super-pixel  vectors  is: 


and  the  breakdown  of  the  super-pixel  zm  into  L  hyper-pixels  is: 

_  /  ~  T  ~  T  \T 

zm  —  \Z{m-l)L+l'>  Z(m-l)L+2->  *  *  *  ?  ZmL)  ? 

(Zq((m—  1)L+1)  ?  Zq((m— l)L+2)?  *  *  *  ■?  ^(mL)) 
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Similarly,  the  breakdown  of  the  conditional  mean  /j,^x  into  M  super-pixel  components  is: 

_  (  T  T  T  \T 

t^z\x  \^Zi \xt  I^Z2\xl  •  *  *  5  I^Zm\x)  ’ 

where  the  rrfi1  conditional  mean  super-pixel  (derived  from  Equation  (4.26))  is  given  by: 


^Zm\x  ~  P'Zm  +  BmSl 


SLBmS[  +  a2IL 


i -i 


(xm  ) , 


(4.68) 


which  may  be  broken  down  further  to  hyper-pixels  resulting  in  Equation  (4.48).  The  quan¬ 
tities  appearing  in  Equation  (4.68)  are  defined  by: 


SL 

=  h®  sT, 

X 

■,xm)T 

P'z 

II 

s:  tT 

■ft 

S' 4 

c, 

©1 

II 

SL  is  L  x  LP,  S  =  SQT  =  IM 


(Vzm  is  LP  x  l)  , 

'  L 

Bm  =  is  LP  x  LP  |  . 

3= 1 


Si 


The  inverted  quantity  in  Equation  (4.68)  reduces  to: 

L 

SLBmSl +  apL  =  0(srb>„,./s  +  ^). 

3  =  1 

Therefore,  the  computation  of  Hzm\x  reduces  to  inverting  diagonal  matrices. 

The  computation  of  zm  using  Equations  (4.64),  (4.67),  and  (4.68)  may  be  accomplished 
in  parallel  by  assigning  a  different  set  of  m  to  each  processing  unit.  This  means  that  no 
parallel  matrix  algorithms  need  be  considered,  making  coding  a  task  considerably  simpler 
than  it  would  otherwise  be.  It  is  important  to  remark  that  Equations  (4.64)  and  (4.67) 
apply  even  if  the  relationship  x  =  Sz  +  rj  is  not  assumed.  This  facilitates  the  construction 
of  a  single  computer  code  that  applies  to  the  case  where  x  =  Sz  +  rj  is  assumed  as  well  as 
to  the  case  where  it  is  not. 


We  now  assume  the  special  case  where  the  P  x  P  diagonal  blocks  of  C%  (and  hence  C~\x) 
are  identical  within  super-pixels.  Thus,  for  each  m  =  1,2 , ,M: 


L  m .  1  Bm2 

Gm,  1  Lrm2 


Bm,L  i 
Gm,L- 


Under  the  special  case  condition,  we  may  characterize  the  definiteness  of  the  inverted  matrix 
of  Equation  (4.58)  in  terms  of  an  and  a,r  For  this  purpose,  we  establish  the  following  lemma. 


Lemma  4.11.1  Given  any  nonzero  vector  v  of  length  P,  the  matrix  T 
positive  definite  if  e  >  0,  and  nonnegative  definite  if  e  =  0. 


T 

VV 

dp  -  T  -  IS 

v1  v  +  e 
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Proof. 


First  assume  that  e  >  0  and  let  T 


11*11* 


vv 


vTv  +  e 


T 

VV 

vTv  +  e 

II  Til 

\\vv  ||f 
vTv  +  e 


so  that 

T 

V  V 

vTv  +  e 


<  1. 


where  1 1|  p  denotes  the  Frobenius  norm.  It  follows  that  T  =  Ip  —  is  nonsingular  and  the 
inverse  of  T  is  given  by  the  Neumann  series: 

OO 

r-1  =  vb", 

n= 0 


For  any  nonzero  vector  x  of  length  P, 


xTtyx 


T  T 
X  VV  X 


V1  V  +  6 


T 

X  V  2 

vTv  +  e 


>  0, 


which  proves  T  is  positive  definite.  Since  4/  is  symmetric,  it  follows  that  T"  is  positive 
dehnite  for  all  nonnegative  integers  n.  which  proves  from  the  Neumann  series  that  T_1  and 
hence  T  is  positive  definite. 


Since  a  matrix  is  positive  definite  if  and  only  if  its  eigenvalues  are  positive,  and  the 
eigenvalues  of  a  matrix  are  continuous  functions  of  its  entries,  it  follows  that  the  eigenvalues 
of  r,  are  positive  or  zero  when  e  =  0.  This  proves  that  T  is  nonnegative  dehnite  if  e  =  0. 


The  fact  that  a  zero  eigenvalue  of  T  actually  occurs  when  e  =  0  is  established  by  observing 
that  vTv  is  an  eigenvalue  of  vvT  with  eigenvector  v,  making  det (vTvIp  —  vvT)  =  0.  Therefore 
T  is  singular  and  hence  has  at  least  one  zero  eigenvalue  when  e  =  0.  ■ 


Theorem  4.11.2  Under  the  special  case  condition,  the  (inverted)  matrix  WCz\xWT  +  Cn  of 
Equation  (4.58)  is  positive  definite  if  either  an  >  0  or  av  >  0  and  is  singular  if  an  =  av  =  0. 


Proof.  From  Equation  (4.62),  it  is  equivalent  to  prove  that  Gm  +  rr(lp  is  singular  if 
&n  =  &ri  =  0  and  positive  dehnite  otherwise.  Under  the  special  case  assumption,  Equation 
(4.61)  becomes: 

Gm  =  ntmGm,i  where 


In  Section  4.7,  we  proved  that  Gm, i  is  singular  when  av  =  0.  It  follows  that  Grn  +  a^Ip  is 
singular  when  an  =  =  0.  Therefore,  we  turn  to  the  case  where  either  an  >  0  or  an  >  0. 


From  Equation  (4.53), 
Gm  +  cr(jp 


'Ym.Gm,  1  T 

ddm,  lSS  Bm\ 


r)m 


Bm.l 


sTBm  is  +  ol 


B 


1/2 
m,  1 


7 m  \  Ip  ~ 


VV 


VTV  +  cr^ 


(JnIP, 

~  a2nBm,l 


B 


1/2 
m,l  ? 
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where  v  =  B^\  s.  Since  the  sum  of  a  positive  definite  matrix  with  a  nonnegative  definite 
matrix  is  positive  definite,  Lemma  4.11.1  implies  that  the  matrix  quantity  within  the  square 
brackets  of  the  last  expression  is  positive  definite  when  either  crn  >  0  or  av  >  0.  From  this, 
it  easily  follows  that  Gm  +  rr2nIp  is  positive  definite  and  the  theorem  is  proven.  ■ 

Experience  with  a  computer  implementation  of  the  Form  1  MAP  algorithm  shows  that 
the  matrix  WCz\xWT  +  Cn  (regardless  of  whether  x  =  Sz  + 1 7  is  assumed  and  regardless  of 
the  special  case  condition)  may  be  poorly  conditioned  when  an  =  0.  In  fact,  taking  an  >  0 
is  a  mechanism  one  can  use  for  diagonally  loading  WC;\XWT  +  CJn .  When  diagonal  loading 
is  used  in  this  manner,  better  performance  of  the  algorithm  is  observed,  as  measured  by 
signal  to  noise  ratio  determined  by  comparing  the  true  z  (assuming  it  is  known)  with  the 
Form  1  MAP  estimate  z. 

Although  WCz\xWT  +  Cn  is  singular  under  the  special  case  and  noise  free  assumptions, 
that  singularity  cancels  when  combined  with  CqxkFT  via  Equation  (4.57).  The  remainder 
of  this  section  derives  an  expression  for  the  Form  1  MAP  estimator  under  the  special  case 
and  noise  free  assumptions.  We  will  see  that  no  matrix  inversions  are  required. 


Under  the  noise  free  assumption  (an  =  av  =  0),  we  see  from  Equation  (4.65)  that 


A $  —  —  (0m  ®  Ip)  > 

7m 


(4.69) 


and  Equation  (4.67)  becomes: 


Ilp - ®  Ip)  (<^m  ®  Ip) 

Km 


M Zm\x  rn  ®  Ip)  Un 


Kn 


Applying  the  properties  of  the  Kronecker  tensor  product  listed  in  the  beginning  of  this 
chapter  yields: 

t^Zm\x  d  ®  Ip)  Umi 

7 m 

which  is  applicable  even  if  the  model  x  =  Sz  +  t]  is  not  assumed.  Similarly,  from  Equation 
(4.68)  we  have: 

Vzm\x  =  Vzm  +  1  (h®  Bm>1  s)  (xm  -  (lL  <g>  sT)  n~Zm)  , 

dm 

where  am  =  srUmls.  Substituting  the  latter  expression  into  the  former  and  collecting  terms 
yields  the  super-pixel  breakdown  of  the  Form  1  MAP  estimator  under  the  noise  free  and 
special  case  assumptions: 


r _ Jj_  ->T 

Km 


Zm.  = 


(J-L  ^rn^rn)  ®  (-^P  Bm  iSS  ) 

Km 


1 


Ih 


Km 


m^m 


1 


2m 


Xm  T  (oJm,  ®  Ip)Vr 1 
7 m 


(4.70) 
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Since  this  expression  is  free  of  matrix  inverses,  it  is  a  candidate  for  real  time  imple¬ 
mentation.  It  is  expected  that  it  will  provide  a  better  result  than  interpolation  of  y  alone 
(without  use  of  x)  while  keeping  complexity  to  a  level  suitable  for  real  time  implementation. 
This  expectation  is  based  on  the  assumption  that  estimates  of  the  covariance  matrix  are 
performed  in  parallel  in  an  efficient  manner. 


4.12  Lagrange  Multiplier  Optimization 


In  this  section,  we  present  an  alternate  derivation  of  the  Form  1  MAP  estimate  using  the 
optimization  technique  of  Lagrange  multipliers.  The  derivation  applies  only  to  the  case 
of  no  noise  present  in  the  observed  data.  Rather  than  an  unconstrained  minimization  of 
(4.29),  we  wish  to  minimize 

C(z)  =  \^z  -  Rq x)TC^x(z  -  pAx)  (4.71) 


subject  to  the  constraint  that  y  =  Wz.  Multiplying  out  the  terms  and  keeping  only  those 
that  are  functions  of  z,  we  arrive  at  an  equivalent  cost  function 


C(z) 


(4.72) 


Thus,  we  are  faced  with  a  minimization  of  a  quadratic  form  with  linear  constraints.  We 
can  solve  the  problem  using  the  method  of  Lagrangian  multipliers  [17].  The  Lagrangian  is 
formed  as 

L(z,  A)  =  \zTC~{xz  -  z  +  XT(Wz  -  y ),  (4.73) 

where  A  is  an  MP  x  1  vector  of  Lagrangian  multipliers.  The  gradient  with  respect  to  z  is 
given  by 

VzUz,X)  =  C;[lz-(C;[1J^x  +  WTX.  (4.74) 

Setting  this  equal  to  zero  and  solving  for  z  yields 


2  =  Pz\x  ~  CZ\XWT\.  (4.75) 

Now  let  us  find  the  A  that  allows  our  solution  to  meet  the  linear  constraint.  To  do  so  we 
impose  the  linear  constraint  on  (4.75)  and  solve  for  A.  This  gives  us 


y  =  Wz  =  W  (**s|l  -  CA,WTX)  . 

Solving  for  A  we  get 

A  =  (WC,lrWT)~'  (Wy.Ax-  y)  . 


Plugging  this  into  (4.75)  yields  our  final  solution 

z  =  y.Ax  -  CzlzWT  [(WCAzWTYl  (W'Xl  -  y) 

This  result  matches  Equation  (4.57)  when  the  covariance  of  the  noise  is  set  to  zero. 


(4.76) 

(4.77) 

(4.78) 
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4.13  Generalization  to  Multispectral  Images 


Up  to  this  point,  the  Nx  1  vector  x  has  represented  an  nvxrih  (i.e.  panchromatic)  image.  We 
now  generalize  the  most  important  equations  (those  required  for  computer  implementation) 
to  the  case  where  x  has  additional  spectral  content.  In  particular,  we  assume  that  x  is  an 
nv  x  rih  x  v  data  cube  written  as  an  Nv  x  1  vector.  In  reference  [18],  the  quantity  v  is 
denoted  by  Q.  We  chooses  a  different  variable  name  in  order  to  avoid  confusion  with  the 
NP  x  NP  permutation  matrix  Q. 

We  adopt  the  following  ordering  and  notation  for  the  Nv  x  1  multispectral  vector  x: 

-r  =  IW1)-  r(2)-  .  TM1 


where  for  1  <  k  <  v  is  an  N  x  1  vector  whose  entries  are  denoted  by 

r(k)  =  (Jk)  (k)  (kpT 

Jy  I  5  *^2  ?  *  *  *  ?  ^  N 


Thus,  contains  the  spatial  content  of  x  corresponding  to  multispectral  band  k.  When 
thought  of  as  an  nv  x  rq,  matrix,  the  entries  of  x(k)  are  stored  in  column  order. 

Our  existing  MATLAB  code  implements  the  Form  1  MAP  algorithm  through  the  use 
of  Equation  (4.67).  Fortunately,  the  generalization  to  multispectral  data  cubes  does  not 
alter  this  equation.  However,  under  the  assumption  x  =  Sz  +  rj,  the  matrix  S  and  the 
vector  rj  take  on  new  meaning,  and  this  has  an  impact  on  how  the  various  terms  and  factors 
of  Equation  (4.67)  are  computed.  Also  affected  is  the  least  squares  approach  presented  in 
Section  4.6  for  estimating  the  spectral  response  matrix  from  x  and  y. 

The  matrix  S  changes  from  an  N  x  NP  matrix  to  an  Nv  x  NP  matrix  and  y  changes 
from  aniVxl  vector  to  an  Nv  x  1  vector.  In  accordance  with  S,  the  P  x  1  spectral  response 
vector  s  becomes  a  P  x  v  matrix  whose  entries  are  nonnegative  and  whose  column  sums  are 
unity. 

For  k  =  1,  2, . . . ,  v,  let  s ^  denote  the  k^1  column  of  the  P  x  v  matrix  s.  Define  the 
N  x  NP  matrix  S ^  by: 

N 

S {k)  =  0[s(fc)]T  (4.79) 

71=1 

=  IN  ®  [s(fc)]T.  (4.80) 

Now  define  the  Nv  x  NP  matrix  S  by  vertical  concatenation  of  the  S^: 


S 


g (!)  ;  gf2).  _  _  _  ■  gU) 


(4.81) 


The  Nv  x  NP  matrix  S  for  the  multispectral  case  is  given  by  S  =  SQ. 
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One  consequence  of  changing  s  to  a  P  x  u  matrix  is  that  the  algorithm  for  its  estimation 
as  given  in  Section  4.6  changes  slightly.  We  define  Xm  analogous  to  Equation  (4.40)  as 
follows: 


x 


(k) 


jL  (o 

j= 1 


In  analogy  to  Equation  (4.41) we  set: 

Xm  =  {ym,S{k))  (1  <  m  <  M). 

In  matrix  form  (analogous  to  Equation  (4.42))  this  becomes: 

where  =  {x^\x£\  ■  ■  ■ ,  %m  )  and  ^  is  the  same  as  before.  Thus,  there  are  now  v  systems 
of  equations  to  solve  in  the  least  squares  sense.  Each  system  is  subject  to  the  constraints 
that  s(k'1  contain  nonnegative  entries  that  sum  to  unity.  The  MATLAB  function  LSQLIN 
may  still  be  used,  but  now  it  appears  within  a  loop  over  k. 

Equation  (4.67),  which  we  used  for  programming  the  Form  1  MAP  algorithm,  ultimately 
depends  on  the  P  x  1  conditional  hyper-pixel  means  ji^ \x  and  the  P  x  P  conditional  covari¬ 
ance  matrices  Gmj  (i.e.  the  j^1  P  x  P  block  of  Cz\x  corresponding  to  the  nfi1  super-pixel). 
Under  the  assumption  x  =  Sz+rj,  these  are  estimated  by  their  (unconditional)  counterparts 
/i£t  and  Bm  :)  through  use  of  Equations  (4.48)  and  (4.53).  However,  these  equations  were  de¬ 
rived  under  the  assumption  of  x  being  panchromatic  and  do  not  apply  to  the  multispectral 
case. 

To  derive  the  generalized  versions  of  Equations  (4.48)  and  (4.53),  we  start  with  the 
more  fundamental  Equations  (4.26)  and  (4.27).  Let  D  denote  the  Nv  x  Nv  matrix  in 
square  brackets  that  appears  in  both  of  these  equations: 

D  =  SCstF  +  Cr,, 

where  Cv  =  cA  INv.  Our  first  problem  will  be  to  invert  D. 

From  Equation  (4.81),  D  is  naturally  partitioned  into  u2  N  x  N  blocks,  where  block 
(i,j)  is  defined  by: 

Ditj  =  S^C~Z[S^]T  +  5iJo3vIN, 

and  Sij  denotes  the  usual  Kronecker  delta  function  (1  if  i  —  j  and  0  otherwise). 

Equations  (4.80)  and  (4.46)  yield  the  following  expression  for  Dt  J: 

N 

D't.j  =  (J)  Ai(b  j)  i 

n=l 

dn(i,j )  =  [s{t)]T Bns{j)  +  Sij  o%. 
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It  follows  that  for  each  n  =  1.  2, ....  Ah  the  u  x  v  matrix  dn  is  given  by: 

dn  =  sTBns  +  <7^  Iu.  (4.82) 

Note  that  these  matrices  reduce  to  the  scalars  di  of  Section  4.7  for  the  case  of  v  =  1. 

Denote  the  Nv  x  Nv  inverse  of  D  l>y  U  and  partition  it  in  a  manner  consistent  with 
D,  using  Uij  to  denote  its  (i,j)  N  x  N  block.  We  now  assume  that  the  zero/nonzero 
structure  of  U  is  identical  to  that  of  D.  Our  assumption  will  be  justified  if  the  resulting 
nonzero  entries  of  U  are  uniquely  determined.  Accordingly,  assume  UtJ  is  given  in  terms  of 
unknown  scalars  un(i,j )  by: 


N 

Uij  = 

71=1 

and  compute  the  un(i,j)  such  that  DU  =  In„.  The  matrix  equation  DU  =  l^v  is  equivalent 
to: 


^  Dj,kUk,j  &i,j  Ini 
k= 1 

where  i  and  j  run  over  1,  2, . . . ,  v.  Given  the  assumed  form  of  U^,  this  becomes: 


N 

© 

77=1 


^  ]  dn (f  >  k)un  ( k ,  j ) 


k= 1 


—  $i,j  In- 


This  however  is  equivalent  to 


^  ^  dn  (b  k^)Un  (fc,  j)  di^ji 

k= 1 

as  n  runs  over  1,  2. ....  Ay  which  is  just  another  way  of  asserting  that  the  v  x  v  matrices 
dn  and  un  are  inverses  of  each  other.  Thus,  the  elements  of  un  are  uniquely  determined  by 
un  =  d~x  and  the  inversion  of  D  reduces  to  the  inversion  of  the  uxu  matrices  di,  d2, . . . ,  djv- 
It  is  seen  from  Equation  (4.82)  that  the  condition  of  dn  will  improve  as  a ^  increases.  Just 
as  the  noise  variance  a \  may  be  considered  a  diagonal  load  to  improve  the  condition  of 
WCz\xWT ,  the  noise  variance  may  be  considered  a  diagonal  load  to  improve  the  condition 
of  SC~ZST. 

Continuing  with  the  generalization  of  Equations  (4.48)  and  (4.53),  we  derive  expressions 
in  terms  of  the  scalars  un(i,j),  which  we  now  know  how  to  compute.  We  proceed  by 
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dissecting  the  factor  STD  1S  appearing  in  Equation  (4.27): 

STD~lS  =  J2J2[S{i)]TUijS{j) 

j  =  1  *=1 

v  v  /  N 

=  EE  ©*<v  (i,j)[sb)]T 

j  —  1  2=1  \  71=  1 

22=1  \  j  =  l  2=1 

Substituting  this  factor  into  Equation  (4.27)  yields: 

C-Z\x  =  Cz-Cz^D-^Cz 

=  ©  (sn  -  Bn^Y/un(i,j)s®[s®]TB 

22=1  \  j  =  1  2=1 

N  /  v  v 

=  ©  [Bn  ~  YYU^^B^B^T 

22=1  \  2=1  j  =  1 


th 

It  follows  that  the  conditional  covariance  matrices  Gmj  representing  the  jLli  P  x  P 
th 

matrix  within  the  super-pixel  of  C% \x  is  given  by: 

V  V 

GmJ  =  BmJ-YY(^i^B^{i))[Bmjs^]T,  (4-83) 

2=1  k= 1 


which  generalizes  Equation  (4.53). 

We  conclude  with  deriving  a  generalization  of  Equation  (4.48).  Dissecting  the  factor 
CZST  D~l  of  Equation  (4.26),  we  have: 


C-zStD~ 1 
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Similarly,  the  factor  a;  —  of  Equation  (4.26)  is: 
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(®(1)  -  S(1)/*5;  x(2)  -  S(2)/a5;  •  •  • ;  ^  -  S{u)n^J 
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Putting  the  factors  together  yields: 
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Therefore,  through  the  use  of  Equation  (4.26),  the  n^1  Pxl  conditional  mean  high  resolution 
hyper-pixel  is  given  in  terms  of  the  corresponding  (unconditional)  mean  high  resolution 
hyper-pixel  by: 


»zn\x  =  Vzn+'52'52  (Xn]  ~  [s{k)]T V'zn)  ««(*,  k)BnS(l\  (4.84) 

2=1  k= 1 

which  generalizes  Equation  (4.48).  Equations  (4.82),  (4.83),  and  (4.84)  are  some  of  the  more 
important  equations  relevant  to  programming  a  generalized  MAP  algorithm  that  applies  to 
both  multispectral  and  panchromatic  inputs. 


4.14  Estimation  of  Conditional  Parameters 

We  now  discuss  our  approach  to  estimating  conditional  statistical  parameters  Cz\x  and  fjtz^x 
under  the  assumption  that  the  linear  model  x  =  Sz  +  rj  does  not  hold.  In  this  section  and 
in  Section  4.13,  we  generalize  the  single  band  panchromatic  array  x  to  a  multispectral  array 
with  v  bands.  Thus, 


/  ->T  ->T  ->T  \  T 

x  =  (x1,x2,...,xN)  , 

where  each  xn  is  a  v  x  1  vector  (for  a  panchromatic  image,  v  —  1). 

Estimating  the  full  NP  x  NP  covariance  matrix  CZ\X)  used  in  Equation  (4.28),  is  imprac¬ 
tical  for  a  typical  size  hyperspectral  image.  To  make  the  problem  manageable,  constraints  on 
the  form  of  this  covariance  are  required  to  bring  down  the  number  of  statistical  parameters 
to  be  estimated.  While  there  may  be  numerous  ways  to  accomplish  this,  we  believe  that  a 
reasonable  approach  is  to  model  the  unknown  high-resolution  hyper-pixels  as  conditionally 
independent  spatially,  yielding 


N 

Pr(z\x^  Pr(KZn\xn)’) 

n=l 


(4.85) 
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which  is  the  multispectral  version  of  Level  1  conditional  independence  presented  in  Section 
4.9.  Writing  out  the  individual  conditional  PDFs  yields 


N 


Pr(z\x)  =  J] 


exp  <!  ~^(zn  -  Vzn\Xn)TCz]x(zn  -  vZn\Xn)  !>,  (4-86) 


Zn\Xn '  Zn\Xn 


Zn  |  Xn  ) 


where  / iZn\Xn  =  E{zn\xn}  and  CZn\Xn  is  the  PxP  covariance  matrix  for  zn  given  xn.  The  rela¬ 
tionship  between  the  individual  hyper-pixel  statistical  parameters  and  the  global  statistical 
parameters  in  Equation  (4.28)  is  given  by 
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and 


P'zjx  \_^zi\x\i  ^Z2\x2i  ’  '  '  ^Z]sr\xi^\  '  (4.88) 

Furthermore,  applying  the  results  in  Equations  (4.21)  and  (4.22)  on  the  hyper-pixels  yields 


Pzn\xn  —  ^  }  4“  CZn,xnCXn  [xn  E{xn}]  (4.89) 


and 


CWk  =  -  CZn,XnC-^CTZntXn.  (4.90) 

Note  that  the  covariance  matrix  of  the  joint  random  vector,  i\)n  =  .7^] T ,  is  related  to 

the  cross-covariance  matrices  in  Equations  (4.89)  and  (4.90)  as  follows 


CL  = 


n  nT 

Xn,Xn  ^ Zn,Xn 
CZn,X„  CZnjZn 


(4.91) 


With  the  simplifying  assumption  of  conditional  independence,  Cz\x  is  block  diagonal, 
reducing  the  number  of  statistical  parameters  from  ( NP )2  to  NP2.  That  is,  the  problem 
reduces  to  estimating  the  hyper-pixel  conditional  statistics  in  Equations  (4.89)  and  (4.90). 
The  number  of  statistical  parameters  may  be  further  reduced  by  assuming  some  or  all 
hyper-pixels  as  having  the  same  conditional  covariance.  If,  for  example,  all  hyper-pixels 
have  the  same  covariance,  then  we  have  only  P2  statistical  parameters  in  the  conditional 
covariance  matrix  to  estimate.  An  area  of  future  work  might  be  to  explore  other  forms  for 
the  conditional  covariance  matrix  (implying  other  assumptions  regarding  the  nature  of  the 
high-resolution  hyper-pixels). 

To  estimate  E {zn},  we  propose  using  the  spatially  interpolated  observed  hyperspectral 
imagery,  denoted  ft Zn .  We  have  found  that  spline  interpolation  tends  to  yield  the  best 
results  here.  To  estimate  E{xn},  we  use  a  spatially  smoothed  version  of  the  band  or  bands 
in  x.  We  have  observed  that  the  best  results  are  obtained  when  the  smoothing  is  done  to 
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mimic  the  way  in  which  jxSn  relates  to  zn.  That  is,  we  degrade  xn  with  a  PSF  and  down- 
sampling  factor  similar  to  that  defined  by  W.  This  degraded  image  is  interpolated  using 
spline  interpolation  to  produce  the  estimate  of  E{;cr,},  which  we  will  denote  ft  ,  for  all  n. 
In  this  fashion  it  may  be  said  that  fiXn  relates  to  xn  as  fiZn  relates  to  zn.  These  estimates 
tend  to  track  the  non-stationarity  in  the  mean  exhibited  by  most  natural  images  [21]. 

The  joint  covariance  in  Equation  (4.91)  must  also  be  estimated  in  order  to  get  the 
conditional  statistics  in  (4.89)  and  (4.90).  In  most  cases  it  will  not  be  possible  to  obtain 
statistically  similar  data  at  the  high  resolution  required.  Thus,  we  attempt  to  estimate 
the  required  joint  covariance  from  the  observed  imagery.  To  do  so,  we  will  estimate  a 
joint  covariance  at  the  lower  resolution  of  the  observed  hyperspectral  imagery  and  apply 
it  at  the  higher  resolution.  While  the  joint  statistics  may  differ  at  different  resolutions, 
it  is  hoped  that  there  is  sufficient  symmetry  of  spatial  scale  in  the  statistical  parameters 
to  provide  a  useful  result.  In  particular,  we  artificially  degrade  the  spatial  resolution  and 

[rri  rri  m  -i  J 

x1 ,  x2  , ...,  xM\ 

The  local  means  at  this  resolution,  obtained  using  the  same  method  applied  to  the  original 
resolution,  are  removed  from  x  and  y.  Now,  the  joint  covariance  information  is  estimated. 
One  relatively  simple  approach  seeks  a  single  global  covariance  using  a  sample  covariance 
estimate.  This  covariance  is  used  as  an  estimate  of  C',/,n  for  all  n. 

However,  in  order  to  more  fully  exploit  the  information  in  x,  we  wish  to  capture  the 
changing  joint  covariance  as  the  spectral  content  in  the  scene  varies  spatially.  Since  it  is 
impractical  to  estimate  a  joint  covariance  for  each  hyper-pixel,  we  use  a  simple  clustering 
approach  based  on  vector  quantization.  To  begin,  we  form  joint  vector  =  \xTm  .  yjn]  G 
E u+p  for  m  =  1, 2, . . . ,  M,  where  E"+p  represents  the  u  +  P  dimensional  real  space.  These 
vectors  are  grouped  into  K  clusters  (or  classes)  using  the  Linde-Buzo-Gray  (LBG)  algorithm 
[22],  The  cluster  centroids  define  the  Voronoi  partitions  of  the  spectral  space.  That  is,  a 
given  vector  in  1Z,,+P  space  is  assigned  to  class  k  if  it  lies  closest,  in  a  Euclidean  sense,  to  the 
centroid  of  cluster  k.  For  each  cluster,  the  sample  joint  covariance  is  computed  using  tpm  for 
m  G  flk,  where  is  the  set  of  all  indices  corresponding  to  Voronoi  partition  k.  To  estimate 
the  joint  covariance,  C.^,n ,  we  assign  =  [xp ,  pT  ]  to  a  partition  and  let  the  covariance 
for  this  high-resolution  spatial  position  n  be  the  corresponding  cluster  covariance.  With 
these  covariance  estimates  in  hand,  the  final  MAP  estimate  can  be  formed.  In  some  cases 
we  have  observed  improvement  in  performance  if  we  form  an  estimate  of  the  conditional 
mean  in  Equation  (4.89),  p>Zn\Xn,  and  then  re-classify  the  image  using  r/y(  =  [.xp,  A!,  JT  for 
the  purpose  of  assigning  cluster  covariances  to  each  high-resolution  position. 
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Chapter  5 

MAP  Algorithm  Performance 


Focusing  on  the  use  of  high-resolution  panchromatic  (single  band)  data  to  enhance  hyper- 
spectral  imagery,  we  show  that  the  MAP  estimator  produces  an  estimate  with  enhanced 
sub-pixel  spectral  content  that  is  evident  not  only  in  the  first  principal  component  image, 
but  in  lower  components  as  well.  The  relationship  of  the  proposed  method  to  some  prior 
methods  is  discussed  in  Section  5.1.  Experimental  results  are  presented  in  Section  5.2. 
Finally,  conclusions  are  given  in  Section  5.3. 


5.1  Relationship  to  Other  Approaches 

It  is  interesting  to  explore  the  relationship  between  the  proposed  MAP  estimate  and  some 
previously  proposed  approaches.  Consider  that  if  one  neglects  the  observation  model  for  y 
in  the  MAP  formulation,  the  resulting  cost  function  would  simply  be  the  second  term  in 
Equation  (4.29).  This  would  lead  to  an  estimate  that  is  the  conditional  mean  in  Equation 
(4.89).  This  result  is  essentially  the  same  as  that  derived  by  Nishii  et  al.  [1]  for  Landsat 
Thematic  Mapper  thermal  band  estimation,  given  specific  choices  for  the  estimates  of  E {zn} 
and  E{xn}.  In  particular,  one  must  use  zero-order-hold  (ZOH)  interpolation  (i.e. ,  pixel 
replication)  on  each  band  of  the  observed  hyperspectral  image  to  estimate  E {zn}  and  average 
the  auxiliary  image  pixels  within  the  span  of  each  low-resolution  hyper-pixel  to  estimate 
E{Tn}.  Using  these  estimates  for  E{£„}  and  E  {:?,,}  guarantees  that  the  average  of  the 
estimated  hyper-pixels  within  the  span  of  a  low-resolution  hyper-pixel  will  be  equal  to 
the  low-resolution  hyper-pixel.  Nishii  et  al.  [1]  explore  the  use  of  both  local  and  global 
covariances.  Thus,  in  comparison  to  their  method,  the  proposed  MAP  framework  is  novel 
in  how  it  explicitly  incorporates  an  arbitrary  system  PSF  and  in  how  it  allows  for  various 
statistical  models  and  estimates  of  the  statistical  parameters.  We  will  use  the  method  of 
Nishii  et  al.  [1]  as  one  of  our  performance  benchmarks  for  comparison  in  Section  5.2. 

Another  method  used  as  a  benchmark  is  the  estimator  proposed  by  Price  [2,7].  This 
method  was  designed  to  combine  multispectral  imagery  with  a  panchromatic  auxiliary  image 
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(i.e.,  v  =  1).  For  bands  strongly  correlated  with  the  panchromatic  sensor,  the  estimate  is 
based  on  a  linear  mapping  of  the  panchromatic  image,  yielding 

Zp:n  =  Q>p%l,n  d*  bp ,  (5.1) 

for  n  =  1,  2, TV  and  p  =  1,  2, P.  This  estimate  is  then  scaled  so  that  the  average  of 
the  high-resolution  hyper-pixels  within  the  span  of  a  low-resolution  hyper-pixel  is  equal  to 
the  low-resolution  hyper-pixel.  The  coefficients,  ap  and  bp  are  estimated  with  least-squares 
regression  using  low-resolution  hyperspectral  band  p  and  a  degraded  version  of  the  panchro¬ 
matic  image  (degraded  to  match  the  spatial  resolution  of  the  low-resolution  hyperspectral 
band).  For  weakly  correlated  bands,  a  look-up-table  (LUT)  method  is  employed  [2,7].  The 
LUT  is  created  based  on  the  relationship  between  the  low-resolution  pixel  values  in  a  given 
band  and  the  corresponding  pixel  values  in  the  degraded  panchromatic  image.  Once  the 
LUT  is  generated  it  is  applied  to  the  full-resolution  panchromatic  image  to  form  a  high  res¬ 
olution  estimate  of  the  desired  band.  As  before,  this  estimate  is  scaled  so  that  the  average 
of  the  high-resolution  hyper-pixels  within  the  span  of  a  low-resolution  hyper-pixel  is  equal 
to  the  low-resolution  hyper-pixel. 


5.2  Experimental  Results 

In  this  section,  we  present  a  number  of  experimental  results  in  order  to  demonstrate  the 
efficacy  of  the  proposed  estimator  in  comparison  to  the  benchmark  techniques.  Simulated 
data  are  used  here  to  allow  for  quantitative  performance  analysis.  The  details  of  the  data 
set  are  provided  in  Section  5.2.1.  In  Sections  5.2.2  and  5.2.3,  quantitative  error  analysis  is 
presented  in  the  spectral  band  space  and  in  the  principal  component  space,  respectively. 
Finally,  noise  analysis  is  presented  in  Section  5.2.4. 


5.2.1  Simulated  Data 

The  simulated  data  sets  are  derived  from  a  hyperspectral  image  collected  by  the  Airborne 
Visible-Infrared  Imaging  Spectrometer  (AVIRIS)  sensor  [23] .  AVIRIS  is  a  scanning  disper¬ 
sive  hyperspectral  imaging  sensor  that  flies  on  the  NASA  ER-2  aircraft  at  approximately 
20  km  above  sea  level  with  a  spatial  resolution  of  approximately  6  m  per  pixel.  The  sensor 
collects  224  contiguous  spectral  bands  in  the  range  of  0.4  to  2.5  pm.  The  specific  scene  used 
has  been  collected  over  Yorktown  Virginia  (Flight  F980703T01,  Run  02,  ID  1828000ST23). 

A  256  x  256  portion  of  the  scene  is  used  as  the  true  high  spatial  resolution  hyperspectral 
image  z.  This  is  artificially  degraded  to  form  y.  A  simple  rectangular  detector  model  is 
used  for  the  system  PSF  [24],  In  particular,  the  PSF  is  a  4  x  4  kernel  with  equal  weights 
of  1/16.  The  image  is  subsampled  by  a  factor  of  4  in  both  spatial  dimensions.  This  PSF 
model  leads  to  a  simple  structure  in  the  matrix  W  described  in  Section  4.3.  This  structure, 
combined  with  the  spatial  conditional  independence  assumption,  allows  us  to  process  each 
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(a)  (b) 

Figure  5.1:  Simulated  observed  images  derived  from  AVIRIS  data,  (a)  False  color  image  of 
principal  components  one,  two,  and  three  of  the  low  spatial  resolution  hyperspectral  data, 
(b)  High  spatial  resolution  panchromatic  image. 


low-resolution  hyper-pixel  to  form  a  corresponding  4  x  4  set  of  hyper-pixels  independently 
(after  the  mean  estimates  are  formed  using  interpolation) .  For  imagery  in  the  mid-  and  long¬ 
wave  infrared,  diffraction  effects  tend  to  become  larger  and  can  be  added  to  the  observation 
model  [24] .  The  associated  high  resolution  sensor  in  this  case  is  modeled  as  a  panchromatic 
broadband  imager  {v  =  1).  These  data  are  formed  by  averaging  the  224  AVIRIS  bands  at 
the  original  resolution.  A  false  color  image,  formed  by  mapping  the  first  three  principal 
components  of  the  low-resolution  hyperspectral  data  to  red,  green  and  blue,  respectively,  is 
shown  in  Figure  5.1(a).  The  simulated  broadband  image  is  shown  in  Figure  5.1(b). 

The  eigenvalue  (variance)  associated  with  each  principal  component  of  the  low-resolution 
hyperspectral  data  is  plotted  in  Figure  5.2.  This  clearly  indicates  that  the  vast  majority 
of  signal  power  is  contained  in  the  leading  components.  For  example,  after  20  components, 
the  eigenvalue  has  dropped  by  approximately  five  orders  of  magnitude  from  the  top  com¬ 
ponent.  In  order  to  reduce  the  computational  burden,  we  process  the  imagery  in  the  PCA 
space  in  the  top  twenty  dimensions.  The  lower  204  dimensions  are  processed  using  spline 
interpolation.  The  processed  components  are  then  transformed  back  to  the  original  spectral 
space.  Note  that  due  to  the  nature  of  the  PCA  transformation,  the  estimation  algorithm 
applies  identically  in  the  principal  component  space  as  shown  in  Section  4.10. 


5.2.2  Spectral  Space  Performance  Analysis 

To  quantitatively  assess  the  performance  of  the  MAP  estimator,  we  compare  the  estimates 
with  the  “true”  hyperspectral  image  (the  original  resolution  AVIRIS  image).  Our  image 
fidelity  metric  is  signal-to- noise  ratio  (SNR),  where  “noise”  here  refers  to  estimation  error. 


71 


Figure  5.2:  Eigenvalue  versus  component  number  for  the  low-resolution  hyperspectral  image. 

This  metric  is  computed  as  the  sample  variance  of  the  “desired”  image  divided  by  the  mean 
squared  error  (MSE).  Scaling  the  reciprocal  of  the  MSE  by  the  variance  of  the  desired  image 
is  helpful  in  allowing  one  to  compare  performance  between  bands  with  significantly  different 
signal  powers.  This  is  particularly  useful  in  principal  component  space,  where  power  in  the 
bands  can  vary  by  orders  of  magnitude. 

The  SNR  versus  wavelength  is  shown  in  Figure  5.3  for  the  MAP  estimator  (K  =  16) 1 , 
the  method  of  Nishii  et  al.  (with  global  covariance  statistics)  [1],  Price’s  method  [2],  and 
spline  interpolation.  For  Price’s  method,  both  the  linear  model  and  LUT  approach  are  used 
and  the  best  of  the  two  SNRs  for  each  band  is  reported.  Here  no  noise  is  introduced  to  either 
the  low-resolution  hyperspectral  imagery  or  the  panchromatic  imagery.  The  effects  of  noise 
are  studied  in  Section  5.2.4.  Note  that  significant  improvement  over  spline  interpolation  is 
obtained  in  many  bands  with  all  the  techniques.  Not  surprisingly,  the  bands  with  the  highest 
correlation  with  the  panchromatic  image  tend  to  have  the  highest  SNRs.  The  spectral  band 
estimates  with  very  low  SNR  are  a  result  of  the  original  data  having  very  low  signal  power 
due  to  atmospheric  absorption.  It  can  be  seen  from  Figure  5.3  that  the  MAP  estimate  with 
(. K  =  16)  provides  the  highest  SNRs  for  these  data. 

1K  refers  to  the  number  of  VQ  clusters  as  described  in  Section  4.14 
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Figure  5.3:  SNR  versus  AVIRIS  band  wavelength  for  the  MAP  estimator  (AT  =  16),  the 
method  of  Nishii  et  al.  (with  global  covariance  statistics)  [1],  Price’s  method  [2],  and  spline 
interpolation. 

5.2.3  Principal  Component  Space  Performance  Analysis 

The  spectral  domain  error  analysis  indicates  that  many  spectral  bands  can  be  significantly 
enhanced  with  the  use  of  the  panchromatic  imagery.  However,  it  is  insightful  to  examine  the 
performance  in  the  principal  component  space.  Table  5.1  shows  the  SNR  in  the  first  5  prin¬ 
cipal  components  for  spline  interpolation,  the  method  of  Nishii  et  al.  (with  global  covariance 
statistics)  [1],  Price’s  method  [2],  and  the  MAP  estimator  for  K  =  1  and  K  =  16.  Clearly 
the  top  principal  component  is  dramatically  improved  by  incorporating  information  from 
the  panchromatic  image.  However,  the  lower  components  are  far  more  difficult  to  enhance 
due  to  the  weak  correlation  with  the  broadband  image.  The  MAP  estimator  does  provide 
a  modest  increase  in  SNR  over  spline  interpolation  for  some  of  the  lower  components,  while 
the  benchmark  techniques  have  lower  SNRs  than  that  obtained  with  spline  interpolation. 
Note  also  that  principal  component  substitution  methods  typically  seek  to  enhance  only  the 
principal  component,  and  do  not  enhance  the  lower  components  at  all.  Thus,  we  believe 
that  any  enhancement  in  these  lower  components  is  a  promising  result.  The  use  of  mul- 
tispectral  high-resolution  imagery  (rather  than  panchromatic)  could  provide  the  means  to 
better  improve  these  lower  components.  To  focus  on  these  lower  components,  Figure  5.4 
shows  the  percentage  SNR  improvement  over  spline  interpolation  for  the  various  estimators 


73 


Table  5.1:  SNRs  for  estimates  of  the  top  5  principal  component  images 


Method 

PC  1 

PC  2 

PC  3 

PC  4 

PC  5 

Spline  Interpolation 

5.87 

6.36 

2.86 

4.86 

3.90 

Conditional  Mean 

26.92 

5.69 

2.48 

4.21 

3.62 

(Nishii  et  al.  [1]) 

Linear  Regression  Method 

27.92 

5.81 

2.48 

4.24 

3.36 

(Price  [2]) 

MAP  (K  =  1) 

34.74 

7.42 

2.99 

5.04 

4.44 

MAP  (K  =  16) 

38.97 

8.37 

3.05 

5.24 

4.48 

in  components  2  through  20.  Note  that  the  MAP  estimator  with  K  =  16  generally  shows 
the  most  improvement.  We  believe  that  the  improvement  seen  with  K  —  16  versus  K  =  1 
is  because  more  correlation  is  present  in  the  individual  classes  than  exists  globally. 

False  color  images  formed  with  the  top  three  principal  components  mapped  to  red,  green, 
and  blue  are  shown  in  Figure  5.5.  In  particular,  the  true  high  resolution  hyperspectral 
image  components  are  shown  in  Figure  5.5(a).  Spline  interpolated  components  are  shown 
in  Figure  5.5(b).  The  estimate  using  Price’s  method  is  shown  in  Figure  5.5(c).  Finally,  the 
MAP  estimate  for  K  =  16  is  shown  in  Figure  5.5(d).  An  enlargement  of  the  upper  right 
corner  of  the  image  is  shown  in  Figure  5.6  for  principal  components  two,  three,  and  four. 
We  believe  that  the  MAP  estimates  generally  appear  sharper  than  the  spline  interpolated 
images  and  exhibit  less  prominent  block  artifacts  than  the  estimates  using  Price’s  method 
(an  observation  consistent  with  the  quantitative  analysis). 

The  number  and  locations  of  the  vector  quantization  code  words  that  define  the  Voronoi 
partitions  has  a  significant  impact  on  the  performance  of  the  algorithm.  Figure  5.7  shows 
the  SNR  versus  the  number  of  partitions  for  principal  components  1,  2,  and  3  for  a  larger 
portion  (512  x  1536)  of  the  AVIRIS  scene.  Note  that  for  8  and  more  partitions,  the  SNRs 
are  markedly  higher  than  with  fewer  partitions.  However,  the  SNR  does  not  appear  to 
increase  linearly  with  number  of  partitions,  as  one  might  expect.  First,  the  performance 
of  the  LBG  algorithm  is  somewhat  sensitive  to  the  initial  starting  point.  The  starting 
vectors  are  selected  by  uniformly  subsampling  the  scene.  Thus,  with  different  numbers 
of  partitions,  different  initial  codewords  are  automatically  selected.  Secondly,  the  LBG 
algorithm  is  designed  to  find  the  most  representative  codewords  (spectra).  It  is  not  designed 
to  optimize  the  MAP  algorithm  performance.  However,  choosing  a  partitioning  scheme 
to  maximize  the  MAP  algorithm  performance  is  a  daunting,  possibly  intractable,  task. 
Furthermore,  these  results  will  depend  on  the  specific  scene  used  and  the  relative  quantities 
of  the  spectra  present.  Alternative  methods  for  partitioning  the  observation  space  may  be 
an  interesting  area  for  further  research. 
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Figure  5.4:  Percentage  improvement  in  SNR  over  straight  spline  interpolation  for  estimates 
of  principal  components  two  through  twenty. 

5.2.4  Noise  Analysis 

In  this  section  we  consider  how  noise  impacts  the  performance  of  the  MAP  estimator.  First, 
we  consider  the  impact  of  noise  in  the  observed  hyperspectral  imagery,  and  then,  we  consider 
noise  in  the  panchromatic  imagery.  The  SNRs  of  the  estimates  of  principal  component  two 
as  a  function  of  the  average  SNR  of  the  observed  low-resolution  hyperspectral  bands  are 
shown  in  Figure  5.8.  Here  no  noise  in  the  panchromatic  image  is  introduced.  One  curve 
shows  the  SNRs  for  the  MAP  estimates  with  K  =  16  when  zero  noise  variance  is  assumed. 
Another  curve  shows  the  SNRs  when  the  correct  noise  variance  is  known  and  used. 

The  SNRs  of  the  estimates  of  principal  component  two  as  a  function  of  the  panchromatic 
image  SNR  are  shown  in  Figure  5.9.  Here  no  noise  in  the  hyperspectral  image  is  introduced. 
Note  that  the  spline  interpolator  does  not  depend  on  x,  and  thus,  is  not  affected  by  noise  in 
the  panchromatic  image.  Also  note  that  when  the  SNR  of  the  panchromatic  image  is  low, 
little  improvement  is  possible  as  correlation  with  the  hyperspectral  bands  is  reduced.  With 
very  low  SNR  panchromatic  imagery,  the  correlation  at  the  lower  resolution  is  not  indicative 
of  the  correlation  at  the  higher  resolution.  This  is  because  the  noise  at  the  lower  resolution 
is  reduced  as  the  panchromatic  image  is  artificially  degraded  to  the  lower  resolution.  The 
use  of  pre-filters  for  noise  reduction  could  help  to  mitigate  this  effect. 
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5.2.5  Multispectral  Auxiliary  Sensor 

We  examined  the  performance  of  the  MAP  estimator  using  a  multispectral  auxiliary  sensor. 
Figure  5.10  shows  the  SNR  for  the  MAP  estimator  (K  =  16)  for  a  different  256  x  256  portion 
of  the  AVIRIS  scene.  We  simulated  a  panchromatic  sensor  (averaging  all  AVIRIS  bands), 
a  visible/near  infrared  dual-band  sensor  (averaging  the  visible  bands,  and  the  remaining 
bands) ,  and  finally  a  six-band  multispectral  sensor  with  bands  similar  to  those  of  the  Landsat 
Thematic  Mapper  sensor  (excluding  the  thermal  band).  Note  that  significant  improvement 
is  possible  in  some  of  the  lower  components,  but  the  number  of  significantly  improved 
principal  components  is  approximately  equal  to  the  number  of  bands  in  the  auxiliary  sensor 
in  these  results.  Modest  improvement  in  many  components  is  possible,  however,  as  shown 
with  the  six-band  auxiliary  sensor  results. 

Figures  5.11  and  5.12  show  some  image  results  comparing  panchromatic  enhancement 
with  multispectral  enhancement.  In  particular,  Figure  5.11(a)  shows  the  raw  low-resolution 
principal  components  1,  2  and  3  in  a  false  color  composite.  Figure  5.11(b)  shows  the  same 
same  components  for  the  spline  interpolated  imagery.  The  high  resolution  components  are 
shown  in  Figure  5.11(c).  The  MAP  estimate  with  K  =  16  and  panchromatic  sharpening  is 
shown  in  Figure  5.11(d).  The  dual-band  sharpened  result  is  shown  in  Figure  5.11(e)  and 
the  six-band  sharpened  result  is  shown  in  Figure  5.11(f).  The  images  are  dominated  by  the 
principal  component  so  that  all  look  quite  good.  However,  some  improvement  is  seen  with 
the  dual-band  and  multispectral  auxiliary  sensor  over  the  panchromatic  sensor.  Figure  5.12 
shows  the  same  set  of  results  as  Figure  5.11,  except  principal  components  2,  3  and  4  are 
used  in  the  false  color  composites.  Here  a  more  dramatic  improvement  can  be  seen  when 
the  dual-band  and  multispectral  auxiliary  sensor  are  used. 


5.2.6  Matched  Filter  Analysis 

To  investigate  the  potential  benefit  of  the  processed  imagery  on  spectral  detection  and 
land-cover  classification,  several  matched  filter  experiments  have  been  conducted.  To  begin, 
an  area  within  the  lake  region  of  the  AVIRIS  imagery  has  been  selected  and  the  spectra 
in  this  area  are  averaged  to  form  a  target  spectrum.  A  spectral  matched  filter  using  this 
target  spectrum  is  applied  to  the  high-resolution  imagery,  yielding  the  result  shown  in 
Figure  5.13(a).  To  cast  this  into  a  detection  framework,  this  matched  filter  image  has  been 
thresholded  to  produce  a  truth  mask,  shown  in  Figure  5.13(b).  Using  the  target  spectrum, 
the  matched  filter  result  using  the  MAP  estimate  with  K  =  16  and  panchromatic  auxiliary 
sensor  is  obtained  and  shown  in  Figure  5.13(c).  For  comparison,  the  matched  filter  result 
using  the  spline  interpolated  imagery  is  shown  in  Figure  5.13(d).  It  is  clear  from  the  matched 
filter  results  that  the  matched  filter  using  the  MAP  imagery  does  a  much  better  job  with 
the  boundary  of  the  lake  region  than  the  interpolated  imagery. 

A  receiver  operating  characteristic  (ROC)  curve  for  the  lake  target  is  shown  in  Figure 
5.14.  This  is  based  on  the  truth  mask  shown  in  Figure  5.13(b).  Note  that  the  area  under  the 
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ROC  curve  for  the  MAP  estimate  imagery  is  clearly  higher  than  that  for  the  interpolated 
imagery.  This  demonstrates  a  definite  improvement  in  matched  filter  detection  performance 
using  the  MAP  imagery  for  this  land-cover  class.  The  most  benefit  has  been  observed  for 
large  landcover  classes  that  are  well  represented  in  the  vector  quantization  codebook  and 
the  corresponding  class  statistics.  Rare  spectra  objects  comprising  only  a  small  number  of 
hyper-pixels  in  the  scene  may  not  be  enhanced  by  this  approach,  relative  to  spline  interpo¬ 
lation.  One  possible  approach  to  enhancing  rare  spectra  objects  is  to  hand  select  some  of 
the  codewords  that  define  the  partitioning  and  providing  sufficient  training  data  for  these 
rare  spectra  objects.  If  the  objects  are  only  present  in  a  few  hyper-pixels  in  the  current  im¬ 
agery,  additional  imagery  may  be  required  showing  the  rare  spectra  objects  amongst  various 
backgrounds  likely  to  be  encountered  in  the  real  data  (the  mixing  aspect  is  as  important  as 
the  rare  spectra  itself). 


5.3  Performance  Summary  and  Conclusions 

This  research  effort  led  to  the  development  of  a  MAP  estimation  framework  for  estimat¬ 
ing  an  enhanced  resolution  image  using  co-registered  high-resolution  imagery  from  another 
sensor.  Here  we  have  focused  on  the  enhancement  of  a  hyperspectral  image  using  high- 
resolution  panchromatic  data.  However,  the  estimation  framework  developed  allows  for 
any  number  of  spectral  bands  in  the  primary  and  auxiliary  image.  We  believe  that  the 
proposed  technique  is  suitable  for  applications  where  some  correlation  exists  between  the 
auxiliary  image  and  the  image  being  enhanced.  The  results  with  AVIRIS  imagery  indicate 
that  a  number  of  methods  do  well  enhancing  the  top  principal  component  image,  where 
strong  global  correlation  exists  with  the  panchromatic  image.  The  lower  component  im¬ 
ages  are  much  more  difficult  to  enhance.  Notwithstanding  this,  we  have  demonstrated  that 
the  proposed  estimator  is  capable  of  providing  modest  improvement  in  some  of  these  lower 
components  (something  not  seen  with  the  benchmark  techniques).  The  spatially  varying 
statistical  model  (i.e.,  K  =  16),  using  vector  quantization,  does  provide  some  additional 
performance  gain  over  global  statistics  (i.e.,  K  =  1).  We  believe  this  is  a  result  of  exploiting 
the  higher  correlations  present  within  the  spectral  classes  (correlations  that  get  “washed 
out”  in  the  global  statistics). 

We  believe  that  one  of  the  merits  of  the  proposed  estimation  framework  is  that  it  allows 
for  an  arbitrary  linear  observation  model.  Furthermore,  the  estimation  framework  opens  up 
opportunities  to  improve  upon  these  results  with  the  use  of  more  sophisticated  statistical 
models  and  methods  for  estimating  the  statistical  parameters  for  those  models.  For  example, 
improved  performance  may  be  possible  using  a  spatial-spectral  model  for  the  desired  high- 
resolution  image  (i.e.,  not  assuming  the  desired  hyper-pixels  are  conditionally  independent, 
given  the  auxiliary  image). 

The  big  question  surrounding  this  research  is:  what  benefit  does  this  processing  bring  to 
the  utility  of  the  processed  imagery?  Clearly  the  imagery  is  enhanced  for  subjective  human 
interpretation.  We  believe  that  the  MAP  approach  arguably  represents  the  most  theoretically 
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sound  method  proposed  to  date  for  merging  hyperspectral  imagery  with  an  auxiliary  sensor. 
Thus,  in  any  applications  where  panchromatic  sharpening  is  used  for  human  interpretation, 
our  method  could  offer  improved  performance.  In  addition  to  the  subject  enhancement, 
we  have  demonstrated  that  some  improvement,  if  only  modest,  is  possible  in  several  of  the 
principal  component  images  (not  just  the  leading  component).  This  in  turn,  may  allow 
spectral  based  classifiers  and  detectors  to  yield  superior  results.  We  demonstrated  that, 
with  large  landcover  classes,  boundaries  are  significantly  enhanced.  Admittedly,  rare  spec¬ 
tra  objects  may  not  be  enhanced  (without  modification  to  the  present  algorithm).  In  light 
of  this  observation,  change  detection  may  be  a  particularly  good  application  for  the  MAP 
processed  imagery.  With  enhanced  boundaries  between  landcover  classes,  such  as  the  lake, 
subtle  changes  in  these  boundaries  will  be  detectable  and  more  accurately  quantified.  Fi¬ 
nally,  landcover  classification  and  segmentation  (possibly  using  unsupervised  clustering  for 
example),  may  yield  more  accurate  and  precise  boundaries  between  large  landcover  classes. 
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(a) 


Figure  5.5:  False  color  images  showing  the  top  three  principal  components  for  (a)  the  true 
high-resolution  hyperspectral  image  (b)  spline  interpolated  components  (c)  linear  regression 
method  (Price  [2])  (d)  the  MAP  estimate  with  K  =  16. 
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Figure  5.6:  False  color  images  showing  principal  components  two,  three,  and  four  for  (a) 
the  true  high-resolution  hyperspectral  image  (b)  spline  interpolated  components  (c)  linear 
regression  method  (Price  [2])  (d)  the  MAP  estimate  with  K  =  16. 
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Figure  5.8:  The  SNRs  of  the  estimates  of  principal  component  two  as  a  function  of  the 
average  SNR  of  the  observed  low-resolution  hyperspectral  bands. 


Figure  5.9:  The  SNRs  of  the  estimates  of  principal  component  two  as  a  function  of  the 
panchromatic  image  SNR. 
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Figure  5.10:  SNRs  for  spline  interpolation  and  MAP  estimator  using  simulated  panchro¬ 
matic,  dual-band,  and  six-band  landsat  TM  auxiliary  senors. 
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Figure  5.11:  Principal  components  one,  two,  and  three  for  (a)  low  resolution  hyperspectral 
imagery  (b)  spline  interpolation  (c)  true  high  resolution  imagery  (d)  MAP  estimate  using 
a  panchromatic  auxiliary  sensor  (K  =  16).  (e)  MAP  estimate  using  a  dual  band  auxiliary 
sensor  (f)  MAP  estimate  using  a  six-band  auxiliary  sensor. 
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Figure  5.12:  Principal  components  two,  three,  and  four  for  (a)  low  resolution  hyperspectral 
imagery  (b)  spline  interpolation  (c)  true  high  resolution  imagery  (d)  MAP  estimate  using 
a  panchromatic  auxiliary  sensor  (K= 16).  (e)  MAP  estimate  using  a  dual  band  auxiliary 
sensor  (f)  MAP  estimate  using  a  six-band  auxiliary  sensor. 
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(c)  (d) 


Figure  5.13:  Matched  filter  results  with  target  spectrum  sampled  from  lake,  (a)  Matched 
filter  output  using  the  true  high  resolution  imagery,  (b)  Binary  detection  map  obtained  by 
thresholding  the  matched  filter  output.  Matched  filter  output  using  the  (c)  MAP  estimator 
(K=16)  and  (d)  spline  interpolation. 
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Figure  5.14:  Receiver  operating  characteristic  curve  for  the  matched  filter  using  the  high 
resolution  imagery  (truth),  MAP  estimator,  and  spline  interpolation. 
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Chapter  6 

Mercury  Multi  Computer 
Implementation 


Due  to  the  assumptions  made  in  Sections  4.3  and  4.11  regarding  the  PSF  and  covariance 
matrix  structures,  the  algorithm  for  computing  the  MAP  estimate  z  may  proceed  one  super¬ 
pixel  at  a  time.  On  a  serial  machine,  a  loop  over  super-pixels  is  required  which  is  the  most 
CPU  intensive  portion  of  MAP_Algorithm.m.1  Therefore,  it  is  the  loop  over  super-pixels 
which  warrants  a  parallel  implementation.  Section  6.1  presents  the  speed  improvements 
gained  by  use  of  the  Mercury  system  as  applied  to  the  super-pixel  loop.  An  important 
related  topic  is  the  use  of  single  precision  and  its  effects  on  timing  and  algorithm  stability. 
This  is  discussed  in  Section  6.2  where  diagonal  loading  is  introduced  to  improve  stability. 


6.1  Timing  Study 

All  timing  results  given  in  this  section  correspond  to  the  most  general  (i.e.  slowest)  case 
where  the  P  x  P  diagonal  blocks  of  the  conditional  covariance  matrix  Cs\x  are  different 
from  one  another  (i.e.  vary  within  super-pixels).  Therefore,  as  a  function  of  the  number 
of  PCA  components  chosen,  all  timings  of  the  super-pixel  loop  represent  the  worst  case. 
Due  to  having  discovered  a  new  way  to  speed  up  the  super-pixel  loop  for  the  general  case 
algorithm,  these  timing  results  are  considerably  better  than  what  was  earlier  thought  to  be 
the  case. 

In  our  MAP  implementation,  one  node  was  designated  as  the  controller.  The  controller 
node  receives  input  data  from  the  host  PC,  and  provides  overall  system  control  for  the  other 
nodes.  After  the  other  nodes  receive  partitioned  data  from  the  controller,  all  nodes  then 
perform  the  desired  processing  in  parallel. 


1This  assumes  covariance  inverses  are  required  within  the  super-pixel  loop,  occurring  for  example  when 
the  linear  model  x  =  Sz  +  r]  is  not  assumed  and  VQ  estimation  of  the  conditional  covariance  is  method  2. 


We  implemented  our  application  using  the  Parallel  Acceleration  System  (PAS),  which  is 
a  C  library  provided  by  Mercury.  Residing  as  an  application  layer  over  the  standard  Mercury 
operating  system  (MC/OS),  PAS  provides  a  convenient  mechanism  for  writing  code  that 
scales  to  any  number  of  nodes.  It  also  allows  for  data  to  be  sent  from  the  controller  to 
the  processing  nodes  in  either  a  BLOCKED  or  WHOLE  manner,  with  usee  synchronization 
between  nodes. 

All  inputs  necessary  for  computing  the  MAP  estimate  are  written  to  a  file  by  the  host 
computer  which  is  read  by  the  controller  node.  The  controller  node  then  scatters  the  input 
to  the  worker  nodes  which  (along  with  the  controller)  work  on  their  share  of  super-pixels. 
Since  the  number  of  super-pixels  is  typically  very  large  in  comparison  with  the  number 
of  available  nodes,  super-pixels  can  be  evenly  spread  among  them.  Since  each  super- pixel 
takes  the  same  time  to  process  as  any  other  super-pixel,  the  work  load  will  automatically 
be  balanced. 

Figures  6. 1-6.4  show  timing  results  as  a  function  of  available  computational  nodes.  Since 
the  curves  shown  were  run  on  the  AdapDev  1280,  we  were  able  to  vary  the  number  of 
available  nodes  from  0  to  8.  The  data  plotted  above  0  nodes  is  timing  results  for  the  host 
computer  alone,  which  is  an  Intel  machine  running  at  1266  MHz  under  Windows  2000.  Data 
plotted  above  1  node  corresponds  to  running  the  entire  super-pixel  loop  on  one  node  (the 
controller  node).  Data  plotted  above  8  nodes  corresponds  to  all  8  nodes  computing  the 
super-pixels  MAP  estimates  in  parallel. 


Timing  of  MAP  Algorithm  on  AdapDev  1280 


Figure  6.1:  Timing  of  MAP  Algorithm  on  AdapDev  1280  (100  PCA  Components) 


The  blue  curves  represent  the  total  time  in  minutes  taken  to  run  MAP_Algorithm  as 
determined  by  the  MATLAB  tic  and  toe  functions.  The  green  curves  represent  the  time 
taken  to  perform  the  CPU  intensive  loop  over  super-pixels,  which  is  computed  in  parallel 
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Timing  of  MAP  Algorithm  on  AdapDev  1280 


Figure  6.2:  Timing  of  MAP  Algorithm  on  AdapDev  1280  (110  PCA  Components) 


Timing  of  MAP  Algorithm  on  AdapDev  1280 


Figure  6.3:  Timing  of  MAP  Algorithm  on  AdapDev  1280  (120  PCA  Components) 
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Timing  of  MAP  Algorithm  on  AdapDev  1280 


Figure  6.4:  Timing  of  MAP  Algorithm  on  AdapDev  1280  (130  PCA  Components) 


whenever  the  number  of  nodes  is  two  or  more.  This  does  not  include  overhead  time  taken 
on  the  controller  for  file  I/O  from  and  to  the  host.  However,  it  does  include  other  small 
overhead  times,  such  as  communication  between  the  controller  and  worker  nodes,  and  the 
time  to  initialize  nodes.  The  black  curves  represent  total  file  I/O  time  between  the  host 
and  controller  nodes,  and  the  cyan  curve  represents  the  time  taken  by  the  host  to  estimate 
parameters.  Since  parameter  estimation  is  solely  a  MATLAB  computation  (i.e.  not  paral¬ 
lelized),  and  file  I/O  time  is  strictly  between  the  host  and  the  controller  (i.e  independent 
of  the  number  of  worker  nodes),  we  see  the  cyan  and  black  curves  are  nearly  constant. 
The  blue  curve  is  well  approximated  as  sum  of  the  cyan,  black,  and  green  curves.  That  is, 
the  total  runtime  of  MAP_Algorithm  consists  mainly  of  three  disjoint  processes:  parameter 
estimation,  file  I/O,  and  super-pixel  computation. 

The  red  curves  represent  a  theoretical  lower  bound  on  the  parallelized  super-pixel  loop, 
ignoring  file  I/O  time  on  the  controller  from  and  to  the  host.  The  red  curve  is  computed 
simply  by  dividing  the  green  time  for  a  single  node,  by  the  number  of  available  nodes.  Thus, 
the  green  and  red  data  point  for  1  node  are  identical.  In  each  of  Figures  6. 1-6.4,  the  red 
curve  is  seen  to  follow  the  green  curve  quite  closely,  indicating  that  the  work  load  is  well 
balanced. 

In  Figures  6.1  and  6.3,  we  note  that  the  time  taken  to  run  the  super- pixel  loop  entirely 
on  the  host  is  significantly  longer  than  the  time  taken  to  run  it  on  a  single  compute  node  of 
the  AdapDev  1280.  This  occurs  even  though  the  host  runs  at  1266  MHz  while  each  Mercury 
CE  runs  at  500  MHz.  The  reason  is  the  use  of  single  precision  and  SAL  function  calls  on  the 
Mercury  CE,  which  can  dramatically  speed  up  computations.  However,  in  Figures  6.2  and 
6.4,  we  see  the  reverse:  the  host  runs  faster  than  a  single  CE.  Timing  results  on  Mercury 
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compute  nodes  are  very  sensitive  to  data  lengths  coming  into  the  various  memory  caches. 
Since  the  number  of  PC  A  components  changes  among  the  figures,  so  do  the  data  lengths 
entering  memory  caches.  Therefore,  variation  in  runtime  between  the  host  and  a  single  CE, 
is  almost  certainly  due  to  changing  data  lengths  entering  memory  caches. 

The  green  curves  of  Figures  6. 1-6.4  show  exactly  what  we  would  expect  if  the  work  load 
is  well  balanced  among  the  compute  nodes.  That  is,  the  curves  have  a  1/N  behavior  just  as 
the  theoretical  lower  bound  (red)  curve  has.  Also,  we  see  that  the  red  and  green  curves  are 
always  close,  which  is  further  evidence  of  efficient  load  balancing.  We  may  also  ascertain 
that  work  loads  are  well  balanced  by  examining  Figure  6.5,  where  150  PCA  components 
have  been  processed  on  8  nodes.  In  that  figure,  we  see  a  time  line  of  events  for  each 
computation  node.  In  particular,  the  green  segments  of  the  plot  represent  the  processing 
time  of  a  node,  i.e.  the  time  it  took  to  process  all  of  its  assigned  super-pixels,  not  counting 
any  communication  time  or  idle  time.  The  start  of  a  green  bar  represents  the  time  at  which 
a  node  received  all  data  necessary  to  compute  its  super-pixels,  and  the  end  of  one  represents 
the  time  at  which  the  node  finished  its  last  super-pixel.  Since  the  green  bars  across  each 
node  have  essentially  the  same  length  (about  246  seconds),  we  see  that  the  work  loads  are 
well  balanced. 
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Figure  6.5:  Time  Line  for  150  PCA  components  on  the  AdapDev  1280 

Also  seen  in  Figure  6.5  are  blue  bars  representing  communication  time  and  red  bars 
representing  idle  wait  time.  The  controller  (i.e.  CE2)  spends  about  40  seconds  during  its 
initial  communication  phase.  This  consists  mostly  of  reading  data  from  the  host  computer 
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but  also  includes  a  small  amount  of  time  pushing  data  to  the  local  memory  of  the  worker 
nodes.  Each  of  the  7  worker  nodes  (i.e.  CE3-CE9)  spends  about  the  same  amount  of 
time  waiting  for  the  controller  node  to  finish  its  initial  communication  phase.  After  a 
worker  finishes  processing  its  super-pixels,  it  has  a  short  wait  while  the  controller  pulls  data 
from  its  local  memory.  Since  the  controller  was  the  first  to  finish  its  processing  of  super¬ 
pixels,  it  also  has  a  short  wait  while  the  worker  nodes  to  complete  their  processing.  The 
final  communication  time  of  the  controller  consists  of  a  negligible  time  pulling  data  from 
the  workers,  followed  by  a  much  longer  time  writing  a  file  used  by  the  host  computer.  No 
communication  time  is  seen  at  any  worker  node  because  data  is  pulled  by  the  controller  from 
the  local  memory  of  each  worker,  as  opposed  to  pushed  to  the  controller  by  the  workers. 

In  order  to  generate  Figure  6.5,  a  series  of  time  stamps  are  inserted  into  C  code  for 
running  on  the  Mercury.  These  time  stamps  are  identical  to  those  for  use  with  Mercury’s 
TATLView  post  processing  program.  After  running  the  code,  a  series  of  event  times  are 
written  to  an  ASCII  file.  Although  this  file  can  be  read  by  the  TATLView  program  to  view 
the  resulting  time  line,  we  have  instead  chosen  to  read  and  process  it  with  MATLAB.  The 
result  is  a  plot  such  as  that  given  by  Figure  6.5. 
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Figure  6.6:  Time  Line  for  224  PCA  components  run  on  the  AdapDev  1280 

Figure  6.6  shows  another  example  running  the  AdapDev  1280  with  8  nodes.  By  taking 
the  number  of  PCA  components  to  be  224,  we  have  increased  the  size  of  the  matrices 
to  invert  to  224.  Qualitatively,  the  green,  red,  and  blue  bars  show  the  same  behavior  as 
Figure  6.5.  Quantitatively,  we  see  that  the  super-pixel  loop  took  approximately  445  seconds 
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and  file  I/O  on  the  controller  took  about  132  seconds.  The  entire  run  of  MAP_Algorithm 
including  parameter  estimation  and  file  I/O  on  the  host,  took  about  11.7  minutes.  This 
case  constitutes  the  full  AVIRIS  data  set,  and  is  the  largest  problem  we  have  run  to  date. 

Figure  6.7  shows  the  same  example  (P  =  224,  N  =  2562,  and  M  =  642),  run  on  the  64 
node  Mercury  (Hawk)  located  at  WPAFB.  Each  node  of  Hawk  is  a  G4  Power  PC  running 
at  400  MHz  and  has  128MB  of  memory  per  node.  This  run  was  important  in  order  to  verify 
the  portability  (and  scalability)  of  the  code,  and  to  verify  that  memory  limitations  on  Hawk 
are  not  exceeded.  Memory  is  a  concern  because  the  input  data  for  all  super-pixels  resides 
in  the  local  memory  of  the  controller  prior  to  scattering  it  to  the  worker  nodes.  The  fact 
that  the  program  ran,  is  conformation  that  we  did  not  exceed  the  128MB  upper  limit. 

It  is  seen  in  Figure  6.7  that  each  of  64  nodes  took  approximately  63  seconds  to  complete 
the  parallelized  super- pixel  loop.  From  Figure  6.6,  we  see  it  took  each  of  8  nodes  445  seconds 
for  the  same  computation.  Since  each  machine  processed  642  super-pixels,  the  AdapDev 
1280  processed  approximately  1.15  super-pixels  per  second  per  node,  while  Hawk  processed 
approximately  1.02  super-pixels  per  second  per  node.  The  different  rates  are  explained  by 
the  different  processing  speeds  of  the  nodes:  500  MHz  for  the  AdapDev  1280  and  400  MHz 
for  Hawk.  We  also  see  that  file  I/O  between  the  host  of  Hawk  and  the  controller  CE  is 
considerably  slower  than  that  of  the  AdapDev  1280.  The  host  of  Hawk  runs  under  a  Solaris 
operating  system  with  limited  memory  and  slow  disk  access.  While  Hawk  might  be  useful 
for  testing  a  real  time  implementation  of  a  MAP  algorithm  based  on  Equation  (4.70),  further 
use  of  Hawk  for  the  existing  implementation  is  unwarranted. 


6.2  Effects  of  Single  Precision  and  Diagonal  Loading 

Our  approach  to  developing  code  to  run  on  a  Mercury  system  consists  of  five  basic  steps. 
First  develop  a  working  code  in  MATLAB,  second  write  a  C-version  to  run  under  Windows, 
third  modify  the  Windows  C  code  to  make  (emulated)  SAL  function  calls,  fourth  port  this 
code  to  the  Mercury  for  running  on  a  single  node,  and  fifth  use  PAS  library  calls  to  parallelize 
it  for  use  on  several  nodes.  Each  step  in  the  process  required  checking  the  answers  obtained 
with  the  original  MATLAB  version. 

Of  the  four  Mercury  systems  we  tested,  the  compute  nodes  have  been  G4  Power  PC  7400 
microprocessor  (or  compatible)  with  AltiVec  technology.  In  order  to  take  better  advantage 
of  this  architecture,  running  in  single  precision  is  important.  Our  experiments  indicate  that 
a  5  fold  speed  increase  can  be  expected  in  changing  from  single  precision  to  double  precision. 
There  are  several  reasons  for  this: 

•  Use  of  the  vector  processing  unit  (VPU)  is  restricted  to  single  precision. 

•  AltiVec  language  extensions  are  restricted  to  single  precision. 

•  The  floating  point  unit  is  twice  as  fast  for  single  precision. 
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Figure  6.7:  Time  Line  for  224  PCA  components  run  on  Hawk  using  64  Nodes 

•  Memory  caches  are  twice  as  efficient  (each  element  takes  half  as  much  space). 

•  The  RACEway  is  twice  as  efficient  (half  as  many  bits  per  word). 

Fortunately,  the  use  of  double  precision  is  not  necessary  for  the  portion  of  the  code  that 
was  parallelized.  This  is  true  even  in  fight  of  the  fact  that  poorly  conditioned  matrices  may 
arise  whose  inverses  are  required.  If  nothing  is  done  to  improve  the  conditioning  of  these 
matrices,  the  inverse  computation  becomes  unstable  resulting  in  MAP  estimates  that  do  not 
match  double  precision  MATLAB  code.  However,  the  matrices  to  be  inverted  are  diagonally 
loaded  when  a  positive  value  of  (input  variable  var_y  of  MAP_Algorithm.m)  is  supplied, 
resulting  in  improved  conditioning.  The  single  precision  C  code  running  in  parallel  on  a 
Mercury  and  the  MATLAB  version  on  the  host  PC  yield  nearly  identical  results  in  all  cases 
tested  provided  an  appropriate  diagonal  load  o \  is  applied. 

Applying  a  diagonal  load  does  raise  some  issues.  One  issue  is  that  supplying  a  diagonal 
load  may  conceivably  have  a  negative  effect  on  SNR.  This  fortunately  does  not  occur.  In 
fact,  supplying  a  diagonal  load  has  a  positive  effect  on  average  SNR  as  indicated  in  Figures 
6.8-6.11,  where  all  computations  were  performed  in  double  precision  MATLAB.  Another 
issue  is  how  to  set  the  diagonal  load.  From  the  same  figures  we  see  that  there  are  a  wide 
range  of  values  we  may  use  for  the  diagonal  load,  all  of  which  improve  average  SNR. 

As  shown  in  Figure  6.8  for  example,  we  see  that  any  diagonal  load  greater  than  0  and  less 
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than  or  equal  to  30  results  in  an  improved  average  SNR.  The  optimal  load  as  determined  by 
the  maximum  average  SNR  is  somewhere  near  12.  Although  the  optimal  load  is  not  readily 
determined  a  priori,  our  experimentation  indicates  that  just  about  any  number  between  say 
1  and  50  should  improve  average  SNR  (compared  to  no  loading)  while  providing  acceptable 
conditioning  for  single  precision  matrix  inversion. 


Figure  6.8:  Effect  of  Diagonal  Loading  on  SNR  (25  PCA  Components) 

An  additional  advantage  to  using  single  precision  over  double  precision  is  the  size  of 
the  problem  that  can  be  loaded  onto  the  Mercury  system.  Since  a  single  precision  variable 
takes  half  the  memory  of  a  double,  the  matrices  involved  can  be  twice  as  big.  This  is 
particularly  important  because  the  local  memory  on  a  node  is  limited.  On  the  AdapDev 
1280  and  the  96  node  Mercury  at  WPAFB  (Titan),  it  is  256MB  per  node,  and  on  the  64 
node  system  at  WPAFB  (Hawk)  it  is  128MB  per  node.  Since  the  controller  must  hold  data 
for  all  super-pixels,  the  size  of  the  problem  is  limited  to  the  available  memory  of  one  node. 

In  summary,  the  use  of  single  precision  is  desirable  on  two  fronts:  speed  and  memory. 
Although  its  use  requires  diagonally  loading,  that  requirement  is  desirable  on  its  own  merits 
as  measured  by  improved  average  SNR.  The  only  drawback  to  diagonally  loading  is  the 
additional  user  input  of  specifying  the  load.  However,  our  experimentation  indicates  that 
there  is  a  wide  range  of  acceptable  loads,  any  of  which  will  improve  average  SNR  while 
simultaneously  conditioning  the  matrices  involved  for  a  stable  inverse  in  single  precision. 
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Figure  6.9:  Effect  of  Diagonal  Loading  on  SNR  (35  PCA  Components) 
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Figure  6.10:  Effect  of  Diagonal  Loading  on  SNR  (35  Spectral  Components) 
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Chapter  7 

Stochastic  Mixing  Model:  Interface 
to  the  MAP  Algorithm 


Two  separate  interfaces  are  presented,  one  where  the  linear  model  x  =  Sz  +  rj  is  assumed 
and  one  where  it  is  not.  In  the  former  case,  we  provide  the  expressions  necessary  for 
estimating  the  required  unconditional  means  and  covariances  and  in  the  later  for  the  required 
conditional  means  and  covariances.  In  either  case,  high  resolution  abundance  maps,  end- 
member  mean  vectors,  end-member  covariances,  and  a  PCA  transformation  matrix  are 
required  as  input.  Only  the  case  where  the  linear  model  is  assumed  (discussed  first)  has 
been  incorporated  into  MAP_Algorithm.m. 


7.1  SMM  Interface  Assuming  a  Linear  Model 

The  stochastic  mixing  model  [25-28],  is  similar  to  the  well-developed  linear  mixing  model  in 
that  it  attempts  to  decompose  spectral  data  in  terms  of  a  linear  combination  of  endmember 
spectra.  However,  in  the  case  of  the  stochastic  mixing  model,  the  endmembers  are 
P  x  1  normally-distributed  random  vectors,  parameterized  by  their  mean  vector  m£k  and 
covariance  matrix  C£k  ,  instead  of  deterministic  spectra.  The  endmember  index  k  ranges 
from  one  to  7Ve,  the  number  of  endmembers  assumed  in  the  scene.  A  finite  number,  Nc ,  of 
mixture  classes  are  defined  as  linear  combinations  of  the  endmember  random  vectors.  These 
mixture  class  random  vectors  are  defined  as 

Ne 

=  ^ak(q)ek ,  (7.1) 

k= 1 

where  ak(q)  are  the  endmember  abundances  associated  with  the  mixture  class,  and  q  = 
1,  2, . . . ,  Nc  is  the  mixture  class  index.  The  endmember  abundances  comprising  the  mixture 
classes  are  assumed  to  conform  to  the  physical  constraints  of  being  between  zero  and  one, 
and  summing  to  unity,  and  are  quantized  at  a  specified  level.  This  is  done  by  a  structured, 
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iterative  search  algorithm  that  determines  all  the  possible  combinations  of  quantized  abun¬ 
dances  (i.e. ,  discrete  levels  between  zero  and  one)  for  the  specified  number  of  endmembers 
that  satisfy  the  physical  constraints.  The  result  is  a  finite  and  discrete  set  of  mixture  classes. 

Each  low  resolution  hyper-pixel  ym  of  the  observed  hyperspectral  image  is  then  assumed 
to  be  a  realization  of  a  particular  mixture  class  random  vector.  Thus,  after  the  SMM  algo¬ 
rithm  converges,  ym  is  assigned  to  a  corresponding  mixture  classe  qm ,  for  m  =  1, 2, . . . ,  M. 
From  this  low  resolution  mixture  class  map,  low  resolution  abundance  maps,  a*,(gm )  are 
formed.  That  is,  there  are  Ne  abundance  maps  giving  the  spatial  abundance  distributions 
for  each  endmember.  Treating  ym  as  a  random  vector,  we  have 


Ne 

Vm  ^  (7.2) 

k= 1 


The  endmember  statistics  are  estimated  from  a  subset  of  the  observed  spectra.  Specif¬ 
ically,  the  sample  mean  and  covariance  estimated  are  computed  for  each  endmember  from 
the  observations  ym  for  rn  <S  ilk-,  where  flk  is  the  set  of  spatial  positions  assigned  to  end- 
member  class  k  in  the  SMM  algorithm.  For  more  detail  on  the  estimation  processes  involved 
in  formulating  an  SMM,  see  [25-28]. 

To  get  a  spatially-varying  statistical  model  at  the  high  resolution,  we  begin  by  bilinearly 
interpolating  the  low  resolution  abundance  maps  to  give  the  high  resolution  abundance 
maps  aj;fc,  for  i  =  1,  2, . . . ,  N,  and  k  =  1,  2, . . . ,  Ne.  Next  we  assume  that  the  high  resolution 
hyper-pixels  obey  the  linear  mixing  relationship  using  the  endmembers  estimated  from  the 
low-resolution  data.  Treating  the  high  resolution  hyper-pixel  Zi  as  a  random  vector,  we 
have: 

Ne 

^  ^  ;  (7.3) 

k= 1 


Thus,  the  spatially- varying  mean  and  covariance  for  the  high  resolution  hyper-pixels  are 
given  by: 


mZi 


C2 


Ne 

'y  '  ai,krnek  1 

k= 1 
Ne 

k= 1 


(7.4) 

(7.5) 


To  incorporate  the  SMM  into  the  MAP  estimator  when  the  model  x  =  Sz  +  tj  is 
assumed,  the  means  and  covariances  of  the  pdf  for  Zi  are  given  by  (7.4)  and  (7.5).  Thus, 
what  would  be  needed  from  the  SMM  is  the  endmember  means  and  covariances  and  the  high- 
resolution  abundance  maps  (or  the  low-resolution  abundance  maps,  that  we  could  bilinearly 
interpolate).  The  endmember  mean  vectors  and  covariance  matrices  could  be  provided  in 
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the  spectral  space  or  the  PC  A  space  (or  simply  the  leading  PC  A  subspace).  If  provided  in 
spectral  space,  and  processing  is  to  be  done  the  PCA  space,  they  can  be  transformed  (using 
the  PCA  transformation  matrix). 

If  the  endmember  means  and  covariances  and  abundance  maps  are  provided  in  a  leading 
PCA  subspace,  we  need  the  PCA  transformation  matrix  used  as  well.  This  way  we  are  sure 
to  transform  the  imagery  to  be  processed  in  a  manner  consistent  with  that  used  to  compute 
the  SMM  parameters.  This  is  due  to  possible  variation  in  how  the  Eigen  vectors  of  the 
spectral  data  are  computed  in  performing  the  PCA  transformation.  Estimates  of  the  lower 
PCA  subspace  can  be  done  using  interpolation. 


7.2  SMM  Interface  Without  Assuming  a  Linear  Model 

Incorporating  the  SMM  into  the  MAP  estimator  when  x  =  Sz  +  77  is  not  assumed  requires 
a  somewhat  more  complex  process.  In  this  case,  consider  joint  endmember  random  vectors 
containing  the  hyperspectral  endmember  and  the  corresponding  spectra  for  the  auxiliary 
image.  In  particular,  let  the  joint  endmembers  be  (P  +  v)  x  1  random  vectors 

uk  =  [el,al]T ,  (7.6) 

where  ctk  is  a  v  x  1  random  endmember  vector  associated  with  the  auxiliary  sensor.  The 
joint  endmember  random  vector  has  mean  mUk  and  covariance  CUk  .  The  joint  endmember 
statistics  are  estimated  using  the  sample  mean  and  covariance  from  the  observations  wrn  = 
[y^l,x'^n]T  for  m  £  where  the  xm  are  artificially  degraded  spectra  from  the  auxiliary 
sensor.  That  is,  the  high  resolution  multispectral  sensor  data  are  degraded  to  a  spatial 
resolution  to  match  the  observed  hyperspectral  data.  Now  define  a  joint  random  vector  at 
the  high  resolution  as  follows 


Wi  =  [z[,xf]T.  (7.7) 

We  will  model  this  as  a  linear  combination  of  the  joint  endmember  random  vectors.  As 
before,  the  relative  abundances  will  be  given  from  the  high  resolution  abundance  maps,  ai>k, 
for  i  =  1,  2, . . . ,  N  and  k  —  1,2, ,  Ne.  Thus, 

Ne 

i  ^  '  @ji,k'U'k 

k= 1 

Note  that  we  are  using  joint  endmember  statistics  estimated  at  the  lower  resolution  for 
the  higher  resolution  estimation  problem.  We  are  also  relying  on  bilinear  interpolation 
to  generate  high  resolution  abundance  maps.  Given  the  linear  relationship  in  (7.8),  the 
spatially- varying  mean  and  covariance  for  Wi  is  related  to  the  joint  endmember  class  statis- 


(7.8) 
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tics  by 


Ne 


^  ^  Q'i,k'muic  i 
k= 1 
Ne 

Cw,  =  y ']  o^cUk , 


(7.9) 

(7.10) 


k= 1 


where  i  —  1,2,. 
be  of  the  form 


,  Ah  Because  Wi  is  a  joint  random  vector,  the  means  and  covariances  will 


mw ,  = 


cv 


^Zi,Zi  ^ Xi,Zi 

CZi,Xi  CXi,Xi 


(7.11) 

(7.12) 


The  conditional  mean  vector  and  covariance  matrix  for  each  spatial  position  are  related 
to  the  joint  statistics  extracted  from  the  right-hand-side  of  (7.11)  and  (7.12)  using 


mZi\Xi  =  mZi  +  C^XiC^  \x. 
C  Zi\Xi  =  C  Zi^z 


m. 


Zi,Xi^Xi,Xi  \^l  •"'Xi 

-1 


^  Zi  ,Xi  ^  Xi  ,Xi  ^  Xi  ,Zi 


(7.13) 

(7.14) 


Thus,  to  incorporate  the  SMM  statistics  into  the  MAP  estimator  when  the  linear  model 
relating  x  and  z  is  not  assumed,  we  require  the  joint  endmember  statistics  mUk  and  CUk . 
We  would  also  require  the  high  resolution  abundance  maps.  With  this  information  we 
would  employ  (7.9)  and  (7.10).  Next  we  would  use  the  right  hand  side  of  (7.11)  and  (7.12) 
to  complete  (7.13)  and  (7.14).  Again,  if  the  information  is  provided  in  a  leading  PCA 
subspace,  we  would  require  the  PCA  transformation  matrix  used  to  ensure  consistency. 
Interpolation  would  be  used  for  the  lower  PCA  subspace. 
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Chapter  8 

Misregistration  Effects 


During  the  course  of  hyperspectral  image  collection,  the  high-resolution  broadband  image 
( x )  and  low-resolution  hyperspectral  image  (y)  can  become  misregistered,  which  means  that 
the  spatial  alignment  of  the  two  images  is  no  longer  consistent.  Misregistration  can  occur 
for  reasons  such  as  a  slight  sensor  misalignment  or  noise  from  various  origins.  One  of  the 
primary  assumptions  of  the  enhancement  algorithm  developed  under  this  program  effort 
is  that  the  two  images  are  in  perfect  alignment.  The  purpose,  then,  of  the  examining  the 
performance  of  the  algorithm  under  less-than-ideal  misregistration/misalignment  conditions 
is  discover  the  algorithm’s  robustness  to  this  effect,  given  that  there  exists  a  real  chance  of 
image  misregistration  in  fielded  hyperspectral  sensor  systems. 

There  are  many  types  of  misregistration  that  can  occur  such  as  translation,  rotation, 
stretching,  and  shear,  or  any  combination  thereof.  For  these  misalignment  types,  the  mis¬ 
registration  errors  that  occur  in-plane  are  by  far  the  easiest  to  compensate  for,  because  an 
inverse  transformation  exists  if  the  mapping  is  one-to-one.  Translation  and  rotation  are  in¬ 
plane  transformations  for  which  a  simple  inverse  transformation  usually  exists.  This  inverse 
transformation  can  easily  be  applied  to  undo  the  misregistration,  and  consequently,  bring 
the  images  back  into  alignment.  Conversely,  out-of-plane  transformations  are  very  difficult 
to  deal  with  given  that  the  mapping  is  usually  one-to-many/many-to-one,  and  finding  the 
inverse  transformation  can  be  impossible. 

The  types  of  misregistration  explored  here  are  restricted  to  translation  and  rotation 
(in-plane  transformations),  although  in  fielded  systems  any  type  is  possible.  The  rationale 
for  restricting  the  misregistration  types  is  predicated  on  the  fact  that  if  the  algorithm  is 
sensitive  to  even  these  types  of  errors,  it  will  definitely  be  sensitive  to  other  types  such  as 
shear  or  stretching  given  these  errors  are  generally  more  severe. 

For  the  experiments  presented,  the  enhancement  quality  was  measured  on  a  pixel-by- 
pixel  basis  between  the  original  high-resolution  hyperspectral  image  (z)  and  MAP  enhanced 
hyperspectral  image  (z).  The  error  metrics  were  the  mean  squared  error  (MSE),  the  mean 
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absolute  error  (MAE),  and  the  signal-to- noise  ratio  (SNR)  given  by  Equation  8.1. 

MSE  4EE  (*,«  -  kn)\ 

2=1  72=1 

MAE=iI)Ekn-4n|,  (8-1) 

2=1  72=1 

SNR  =  Ef=i  Ell  (*,»  -  ^)2/Ef=i  Ell  (*,»  -  4»  -  K)2, 

where  /q  =  ^  Eli  A,n  and  /i*  =  £  Eli(^,n  -  4«)- 

The  data  utilized  was  a  single  AVIRIS  image,  where  the  broadband  image  was  created  by 
averaging  the  hyperspectral  bands  of  the  original  image.  The  low-resolution  hyperspectral 
image  was  created  by  applying  a  4x4  averaging  PSF  to  the  original  hyperspectral  image.  For 
comparison,  the  enhancement  algorithm  was  evaluated  against  a  random  broadband  image 
as  well  as  spline  interpolation  of  the  low-resolution  hyperspectral  image.  The  enhancement 
algorithm  was  run  in  PCA  mode,  keeping  P  =  40  principal  components. 


8.1  Rotation  Experiments 

The  first  set  of  experiments  conducted  involved  rotating  the  high-resolution  broadband 
image  a  small  amount  about  the  center  of  the  image  to  bring  the  imagery  out  of  alignment. 
The  rotation  transformation  is  defined  by  Equation  8.2 


'  i' ' 

cos (9)  sin(0) 

i 

.  f  . 

— sin(0)  cos(0) 

.  j  _ 

where  9  is  the  angle  of  rotation,  [i  j]T  are  the  original  coordinates,  and  [i'  j']T  are  the 
coordinates  of  the  rotated  broadband  image.  After  the  image  was  rotated,  the  image  was 
cropped  back  to  the  original  image  dimensions  and  imputed  in  the  enhancement  algorithm. 
Figure  8.1  shows  the  enhancement  results  as  a  function  of  rotation  angle  (9). 

Figure  8.1  clearly  demonstrates  that  even  a  small  misregistration  error  can  reduce  per¬ 
formance  significantly.  After  about  1°  of  rotation,  the  map  algorithm  performs  at  the  level 
of  the  baseline  techniques. 


8.2  Translation  Experiments 

The  next  experiments  were  conducted  for  the  case  of  translation  (both  vertical  and  hori¬ 
zontal  simultaneously)  of  the  broadband  panchromatic  image  relative  to  the  low  resolution 
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map  algorithm  as  function  of  angle 
interpolation  of  y  alone 
— •-  random  x  into  map  algorithm 


Figure  8.1:  Rotation  misregistration 


hyperspectral  image.  The  simple  equation  for  this  transformation  is  given  by: 


i* 

_ 

a 

+ 

i 

3  \ 

a 

_  3  _ 

(8.3) 


where  a  is  translation  offset  factor.  The  results  are  summarized  in  the  Figure  8.2. 

These  results  reflect  the  rotation  results  showing  that  once  the  broadband  image  is 
perturbed  by  even  1  pixel,  the  map  algorithm  is  no  longer  effective.  Even  a  misregistration 
error  of  a  fraction  of  a  pixel  can  degrade  performance  substantially. 


8.3  Misregistration  Conclusion 

Overall  then,  these  results  demonstrate  the  algorithm  requires  well-aligned  and  registered 
images  to  be  effective.  The  effects  due  to  out-of-plane  misregistration  were  not  tested,  but  it 
is  likely  that  algorithm  sensitivity  would  be  comparable  to  the  in-plane  case.  If  it  is  known 
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Total  SNR  ratio  [  var(z)  /  var(z-zhat)  ]  Total  Mean  Squared  Error 


-0-  map  algorithm  as  function  of  translation 
*-  interpolation  of  y  alone 
random  x  into  map  algorithm 


Figure  8.2:  Translation  misregistration 


that  target  sensor  platform  can  fall  out  of  alignment,  it  is  necessary  to  add  some  front-end 
preprocessing  to  register  the  images. 
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Appendix  A 

Analysis  of  the  Form  2  MAP 
Estimate 


As  we  saw  in  Section  4.5,  the  Form  1  and  Form  2  MAP  estimates  are  identical  and  equal  to 
the  expected  value  of  z\tp  under  the  assumption  x  =  Sz  +  rj.  This  appendix  may  therefore 
be  considered  as  redundant,  however  it  does  present  alternate  derivations  that  might  be 
useful  at  the  coding  level.  It  also  serves  as  a  check  on  the  derivations  presented  in  Section 
4.11.  The  three  conditions  of  Section  4.3  are  assumed  and  the  definitions  of  Bm  j  given  by 
Equations  (4.49)  and  (4.50)  are  adopted. 

Starting  with  the  Form  2  MAP  estimate  presented  in  Section  4.5,  this  section  derives 
an  equivalent  expression  that  requires  fewer  matrix  inversions  (from  4  to  2).  Assuming 
the  special  case  where  the  P  x  P  blocks  of  C~  are  identical  within  super-pixels,  we  then 
prove  that  the  matrix  appearing  as  an  inverse  in  the  new  expression  is  singular  if  and  only 
if  an  =  av  =  0.  A  separate  expression  is  developed  under  the  noise  free  and  special  case 
assumptions  and  in  this  expression  which  is  identical  to  that  of  the  previous  section. 

From  Section  4.5,  the  Form  2  MAP  estimate  is  given  by: 

z  =  [C-1  +  WTC~1W  +  STC~1S] _1  [C~xiiz  +  WTC~ly  +  STC~1x\  (A.l) 

Let  R  =  [C~l  +  WTC~1W]  1  and  use  the  matrix  inversion  lemma  (see  Equation  (4.32))  to 
obtain1  : 

R  =  CZ-  CZWT(WCZWT  +  Cn)~lWCz.  (A. 2) 

Apply  the  matrix  inversion  lemma  again  to  obtain: 

[R-1  +  STC~1S ] _1  =  R  -  RSt(SRSt  +  C^SR. 

Comparing  with  Equation  (4.27),  we  see  that  the  matrix  R  is  really  Cz\y. 
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Then 


z  =  [R  -  RSt(SRSt  4  C^SR]  [C;'fiz  4  WTC~1y  4  STC~1x] 

=  [R  -  RSt(SRSt  4  C^SR]  (C7V2  +  WTC~ly) 

4  [R  -  RSt(SRSt  4  C^SR]  ( STC~1x ) 

=  u  4  v  —  RSt(SRSt  4  C'?7)_15'u,  (A. 3) 

where 

u  =  Ft  +  WTC~1y) , 

V  =  [fl  -  RSt(SRSt  +  Cy-’SR]  (STC-1x) . 

After  some  algebra,  new  expressions  for  u  and  v  may  be  derived  2: 

u  =  yz  +  CzWT(WCzWT  +  Cn)-\y-Wnz)  (A.4) 

v  =  RST(SRST  +  Cr])~1x.  (A. 5) 

Substituting  (A. 5)  into  Equation  (A. 3)  yields: 

z  =  u  +  RSt(SRSt  +  Cv)-\x-  Su).  (A. 6) 

Equations  (A. 2),  (A.4),  and  (A. 6)  taken  together  specify  the  Form  2  MAP  estimate.  It 
is  now  straightforward  to  derive  analogous  results  in  terms  of  the  simpler  quantities  W  = 
WQT,  S  =  SQT,  R  =  QRQT,  Cz  =  QCZQT ,  u  =  Qu,  y~z  =  Qyz,  z  =  Qz,  and  *  = 
(xi,x2, . .  .,xn)t: 

R  =  CZ-CZWT(WCZWT +  Cn)~1WCz,  (A. 7) 

u  =  ^  +  CjWtJWCzWt  +  Cn)~\y  -W^z),  (A. 8) 

z  =  u  +  RST(SRST  +  Cv)~\x-Sn).  (A.9) 

Combining  Equations  (A. 7),  (A. 8),  and  (A.9)  and  collecting  terms  yields: 

5  =  (INp-A^W)(INp-eS)y,z  +  (INp-eS)A^y  +  ex,  (A. 10) 

where  the  NP  x  MP  matrix  A{2>  and  the  NP  x  N  matrix  0  are  given  by: 

=  CZWT(WCZWT  4  Cn)-\ 

©  =  RST(SRST  P  CV)  '. 

Comparing  Equations  (A.l)  and  (A.  10),  we  see  that  the  number  of  matrix  inversions 
has  been  reduced  from  4  to  2.  Since  the  two  eliminated  are  inverses  of  Cn  and  Cv ,  this  is 
advantageous  when  either  Cn  or  Cv  is  not  taken  as  a  multiple  of  the  identity  matrix.  In 

2Comparing  with  Equation  (4.26),  we  see  that  the  vector  u  is  really  Hz\v- 
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analogy  with  the  Form  1  MAP  estimate,  the  Form  2  MAP  estimate  first  computes  z  using 
Equation  (A.  10)  and  then  computes  z  from  z  —  QTz. 

Analogous  to  Equation  (4.61)  of  Section  4.11,  let 

L 

Bm  =  ^  (A.  11) 

3= 1 

where  Bm:1l  is  the  P  x  P  matrix  given  by  Equations  (4.49)  and  (4.50). 

The  results  of  Section  4.11  (with  Cz\x  replaced  by  C-,)  prove  that: 

M 

AW  =  CzWT(WCzWT  +  Cn)~l  =  ®A^,  (A. 12) 

m=  1 


where  A 


(2) 


is  the  LP  x  P  matrix  given  by: 


am 


(-^m  “1“  ^77, -^p) 

"F  ^n^p) 


From  (4.10),  (A. 7),  and  (A. 12): 


R  = 


M 

InP  -  ®  AMWm 

m= 1 


cs. 


(A-13) 


(A-14) 


It  is  readily  seen  from  Equation  (A.  14)  that  R  is  block  diagonal  with  LP  x  LP  blocks. 
However,  the  LP  x  LP  matrix  Am  Wm  is  not  block  diagonal  with  P  x  P  blocks  which 
implies  that  R  is  not  block  diagonal  with  P  x  P  blocks.  Nevertheless,  we  will  still  be  able 
to  establish  a  super-pixel  breakdown  of  the  Form  2  MAP  estimate,  valid  under  the  special 
case  assumption  which  will  be  similar  to  the  Form  1  breakdown  given  by  Equation  (4.67). 

Notice  that  Bm  is  positive  definite  since  each  Bm  ]  is  positive  definite.  This  implies  that 
the  matrices  Bm  +  crnIp  appearing  in  Equation  (A.  13)  are  nonsingular  regardless  of  whether 
cr^  =  0  or  an  >  0,  proving  that  WCpVT  +  Cn  is  nonsingular  regardless  of  whether  <x„  =  0 
or  an  >  0.  However,  the  inverse  of  SRS T  +  Cv  also  appears  in  the  Form  2  MAP  estimate  as 
specified  by  Equation  (A.  10).  We  will  prove  below  that  under  the  special  case  assumption, 
the  singular  or  nonsingular  status  of  SRST  +  Cri  depends  on  the  zero  or  nonzero  status  of 
both  cq,  and  an.  With  this  in  mind  we  prove  the  following  preliminary  results,  which  will 
also  be  utilized  when  the  Form  3  MAP  estimate  is  analyzed.  As  before,  jrn  denotes  the 
inner  product  jm. 
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Lemma  A. 0.1  ||Am||i?  <  1  whenever  av  >  0  or  an  >  0  where 

*  _  S  Bm  \  (qVre-®rra,l  P  ffijLp)  1®  /  ->  ->T\ 

Am.  A  y^m^m )  ? 


sTBmjls  +  o* 


and 


denotes  the  Frobenius  norm. 


(A.15) 


Proof.  Since  ym  = 


IIA. 


jT 

rr 


m\\F 


s  Bm  i  (. Bm  i  +  5m/p)  Bm  i s 
sTBmAs  + 


where  Sm  =  cr^/7m.  Adding  and  subtracting  SmIp  gives  the  following: 

Bm, i  (-^m, i  A-  dmIp)  [  dmIp  ( Bm  \  -\-  SmIpf  (Bni  i  T  SmIp) 

=  ~ Are  +  Sm Ip)  1  +  Ip. 

Therefore, 


II A 


m  F  - 


S  Are  T  An-^p)  T  Ip  Bm,  1® 

ST5m,iS  +  0-2 

An®  {Bm,  1  T  dmIp)  Bm  \S  T  S  Bm  \ s 

ST5m>iS  +  0-2 


and  ||Am||ir  <  1  if  and  only  if: 


Are®  (  Bm  ]  T  5mIp )  Bm  \ S  <C  (T^. 

Since  Are  >  0,  we  see  that  this  last  inequality  holds  if  an  >  0.  It  also  holds  if  an  >  0  since 
in  this  case  Sm  >  0.  ■ 


Lemma  A. 0.2  IL  —  Am  is  singular  if  and  only  if  an  —  av  —  0. 

Proof.  By  Lemma  A. 0.1,  it  is  enough  to  prove  that  Ip  —  A.m  is  singular  whenever  an  = 
av  =  0.  Assuming  an  =  av  =  0,  we  have  from  Equation  (A.15)  that 

S  Bm  \  (7m-®m,l)  Bm  \S 


Am  = 


srBm>1s 


(^m)  3 


1 

- L 

Am 


Therefore, 


Am  {.Am^-L  ^m^rn)' 

Am 


Since  ym  is  an  eigenvalue  of  w>ma iff  (with  eigenvector  um),  it  follows  that  the  determinant 
of  Ip  —  Am  is  zero  making  Ip  —  Am  singular.  ■ 

The  following  theorem  establishes  sufficient  conditions  under  which  all  matrices  appear¬ 
ing  as  inverses  in  the  Form  2  MAP  estimate  of  Equation  (A.  10)  are  nonsingular. 
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Theorem  A. 0.3  For  the  special  case  Bm}j  =  Bm: i  (1  <3  <  L,  1  <  m  <  M),  the  matrix 
H  =  SRST  +  Cv  appearing  as  an  inverse  in  the  Form  2  MAP  estimate  is  nonsingular  when 
either  an  >  0  or  a v  >  0,  and  is  singular  if  an  =  av  =  0. 


Proof.  Since  the  matrix  SRST  +  Cv  is  block  diagonal  with  LP  x  LP  blocks,  its  invert- 
ibility  is  determined  from  the  invertibility  of  each  block.  From  Equation  (A.  14),  the  m^1 
diagonal  block  of  SRST  +  Cv  is  given  by: 

Hm  =  SL(ILp-A^Wm)BmS^  +  a2vIL, 

=  SLBrnSl  Pa2IL  -  SLA^WmBmS[,  (A. 16) 


and  we  wish  prove  that  an  =  av  =  0  is  equivalent  to  Hm  being  singular  for  m  =  1,2, ,  M. 
Since  under  the  special  case  assumption, 

L 

R  m  R  w .  /  I L  0  Bm  i, 

3= 1 


it  easily  follows  (using  Sl  =  II  ®  st)  that: 

SLBmS[ 

ShBmSi  +  O'2 1 1 


h  ®  Bm, i  s 
(sTBmAs)  IL, 
(sTBmAs  +  a2)  IL. 


Likewise,  from  Equations  (A.  11)  and  (A.  13): 


Bm  1  with  (jJrnUJrn , 


k  -1 


A$  =  L Bm>1  (7 mBm,l  +  <72nIp) 
it  easily  follows  (using  Wm  —  ^m  0  Ip)  that: 

^m^m  ®  Bm,l  A  CT^/p) 


-1 


Thus, 


^rrt  ^^m-Bm  ®  Bm  \  (7m^m,  1  T  <JnIp)  Sm?x, 

SLA$WmBmSZ  =  c omuTm  0  (sTBmA  (7 mBm,i  +  <? Bm>1  s)  , 

(s  R m .  1  ip/mBm,!  *4"  ®n^p)  -®m,l®) 

and  substitution  into  Equation  (A.  16)  yields: 

Hm  (s  Bm\S  +  crn)  I L  (s  Bm  \  (ftmBm,!  P  ^n4p)  UmUlm. 
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Therefore 


Hm  (s  Brn  i  S  T  (7 (I L  A m )  > 

where  Am  is  given  by  Equation  (A. 15).  By  Lemma  A. 0.2,  the  theorem  is  proven.  ■ 

As  we  will  see  in  Appendix  B,  the  proof  of  Theorem  B.0.5  shows  that  II  —  Am  is  positive 
definite.  This  implies  that  Hm  and  hence  H  =  SRST  +  Cv  is  positive  definite  under  the 
special  case  assumption  and  the  assumption  that  either  an  >  0  or  av  >  0. 

Below,  the  Form  2  MAP  estimate  given  by  Equation  (A.  10)  is  broken  down  into  super¬ 
pixel  components,  making  it  amenable  to  parallel  processing.  By  Theorem  A. 0.3,  it  is  valid 
under  the  special  case  assumption  and  the  assumption  that  either  an  >  0  or  av  >  0.  As  with 
super-pixel  breakdown  of  the  Form  1  MAP  estimate,  the  LP  x  1  super-pixel  components 
are  denoted  by  zm  and  are  given  in  terms  of  the  P  x  1  low  resolution  hyper-pixel  ym .  Also 

appearing  is  the  Lx  1  vector  xm.  The  lack  of  a  tilde  on  xm  is  intentional  since  it  represents 
th 

the  rrr  set  of  L  consecutive  entries  of  x  itself  and  not  a  permuted  version  of  x.  From 
Equation  (A.  10),  the  rrfi1  super-pixel  Form  2  MAP  estimate  is  given  by: 

Zm  =  ( IlP  ~  (iLP  ~  +  (jLP  ~  Vm  +  0m®m( A. 17) 

where 

0  =  RST(SRST  +  C,)-1  =  ®^=1©m. 


From  equation  (A.  14): 

©m  =  {ILp-A^Wm)BmSlH-\ 
=  SL(ILP^A^Wm)BmSj 


Hn 


m 

T  ,  2  T 

l  +  Tpr- 


The  special  case  assumption  provides  significant  simplifications  of  Hm,  0m,  and  A 


(2). 


where 


Hm  = 

(sTBm^S  +  crj)  (IL  ~  A m)  , 

(A.18) 

©m  = 

{Jl  ®  Bm  iS  UJrnUJrn  ®  Cms)  Hm  , 

(A.19) 

4(2)  - 

^ m  ~ 

®  Bm  i  Bm,l  4“  ^7 %Ip) 

(A.20) 

A  m 

S  CmS  / 

sTBm,lS  +  ^  > 

(A.21) 

Cm 

Bm,  1  4“  ^7 %Ip)  Bm,l- 

(A.22) 

(2) 

Alternate  expressions  for  Am  and  Cm  are  more  efficient  computationally: 

A^  =  Ljm  ®  —  [lP  —  5m  (Bm  i  +  5mIp)  x]  , 

'"V  L 

Jm 


Cm  —  (Bm,  1  —  8mIp  +  $m  (Bm,l  +  ^m^p)  )  • 
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where  Sm  =  These  may  be  derived  from  Equations  (A. 20)  and  (A. 22)  by  adding 

and  subtracting  SmIp  from  £>m, i  and  distributing  the  factor  ( Brrt: ,  +  5mIp)~l . 

Equation  (A.  17)  is  valid  whenever  either  an  >  0  or  av  >  0  and  the  P  x  P  blocks  of 
Cz  are  assumed  to  be  identical  within  super-pixels.  As  was  the  case  for  the  Form  1  MAP 
estimate,  the  computation  of  zm  using  Equation  (A.  17)  may  be  accomplished  in  parallel  by 
assigning  a  different  set  of  m  to  each  processing  unit,  thereby  avoiding  the  need  for  parallel 
matrix  algorithms. 

It  was  shown  in  Theorem  A. 0.3  that  Hm  is  singular  when  arj  =  an  =  0,  which  would 
seem  to  indicate  that  the  noise  free  case  is  incompatible  with  the  Form  2  MAP  estimate. 
However,  as  we  prove  below,  the  singularity  represented  in  0m  cancels  in  the  noise  free 
case,  allowing  the  estimate  to  be  expressed  in  a  way  that  no  inverse  matrix  appears.  The 
following  lemma  facilitates  the  argument  used. 


Lemma  A. 0.4  For  any  (3  with  0  <  (3  <  A-, 


1 


(II - ^m)  {Jl  ~  f3flm) 


-1 


I L  ^m-> 

Tm 


where 


Proof.  For  0  <  (3  <  — 

7m 

||/A2m||^  =  1 1 1 1 F  =  <  1, 

proving  that  p  —  (Klm  is  nonsingular.  Therefore,  it  is  equivalent  to  prove  that: 

(P - ^m)  =  (P - ^m)(P  —  (3Clm), 

'Km  dm 

Multiplying  this  out  and  cancelling  terms  shows  that  the  above  condition  is  equivalent  to 
^  being  a  projection  map: 


Qr. 


1 1  r 


um  \  _  u  um 

,  Tro  )  Tm 

which  is  easily  verified.  ■ 

We  will  conclude  the  section  by  deriving  a  single  expression  for  z  that  is  free  of  any  matrix 
inverses  and  applies  under  the  special  case  and  noise  free  assumptions.  To  accomplish  this, 
we  take  an  =  0  in  Equations  (A.  18)  to  (A. 22): 

C  -  —R 


Tm 

q(2) 

CB)  Ip-) 

Tm 

A  m 

dm  (Jl  /dm{ 
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where 


dm  S  Bmi  S  +  (7  , 


7  2 

dm'Hm 


By  Lemma  A. 0.4,  when  crn  =  0  and  av  >  0,  Equation  (A.  19)  becomes: 


*171 

1 


7r? 


®m  7  (  II  ®  ^m^m  ®  ^  lS  )  (-Tl  /3m(lJmCJm)  , 

Qjrn 


d 


m 

1 


1 


dn 


d L  ^m^m,  )  ®  -®ra,lS 

7m 

1 

7m 


(4  -  A 


-»T  \  —  1 

m^m^m  J  5 


d L  ^m^m  )  Pm^m^m) 


1  (  T  1  _T  t  „ 

“7  (  d L  ^m^m  )  ®  -®m,  \  S . 

^m  \  7m 


It  follows  that  as  oy  — >  0: 


9r 


(  II  ^m^m  )  ® 

^  1  f 

Jm 


where  am  —  cT 


ST5miiS. 


Substituting  these  results  into  Equation  (A.  17)  and  simplifying  yields  the  super-pixel 
breakdown  of  the  Form  2  MAP  estimate  valid  under  the  noise  free  special  assumptions: 


(II  -  -A^)  <8>  (Ip  -  —Bmtl ssT) 
7 m 


+ 


1 


7n 


a,  ~ 


B'Zm 

1 


H-  (^m  ®  dp^Um-  (A. 23) 

7 771 


As  a  check,  we  see  that  Equations  (4.70)  and  (A. 23)  are  identical. 
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Appendix  B 

Analysis  of  the  Form  3  MAP 
Estimate 


The  analysis  of  the  Form  3  MAP  estimate  presented  in  this  appendix  is  analogous  to  that 
of  the  Form  1  and  Form  2  MAP  estimate  presented  in  Section  4.11  and  Appendix  A.  We 
first  express  the  estimate  in  terms  of  simpler  matrix  quantities  by  applying  the  permutation 
matrix  Q.  Then,  we  break  down  the  coefficient  matrix  of  the  system  into  super-pixel 
blocks,  and  analyze  the  special  case  where  the  P  x  P  blocks  of  the  covariance  matrix  Cg  are 
identical  within  super-pixels.  The  section  concludes  by  providing  a  simplified  expression  for 
the  Form  3  MAP  estimate  broken  down  into  super-pixel  components,  applicable  under  the 
special  case  nose  free  assumption.  Significant  use  is  made  of  the  properties  of  the  Kronecker 
tensor  product  listed  in  the  Chapter  4. 

Starting  with  Equation  (4.39)  of  Section  4.5,  we  apply  the  permutation  matrix  Q  to 
derive  the  equivalent  expression: 


=  V-z  +  A{3) 


y 

X 


(B-l) 


where 


= 


+ 


-i  -i 


Cn  0 

0  c„ 


Similar  to  the  Form  2  MAP  estimate,  we  examine  the  conditions  under  which: 


H  = 


(B.2) 
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is  singular  or  nonsingular.  The  three  conditions  of  Section  4.3  are  assumed  and  the  defini¬ 
tions  of  Bmj,  and  Bm  given  by  Equations  (4.49),  (4.50),  (A. 11)  are  adopted.  Thus, 


M 


S  —  Im®Sl  — 


m= 1 
L 


SL  =  /l®st  =  ©sT, 

3= 1 


and 


M 


W  -  ©  0  Ip)  . 


m=  1 


From  equation  (B.2), 


H  = 


WCtW^  +  cn  WC~ZST 
rT  sc~zst  +  c, 


sczw 7 


V 


Using 


M 


wczw 


-  Bm, 

m—  1 
M 

SC~ZST  =  0  SLBmSl, 


WCZSJ 


we  break  down  H  into  blocks: 

M 


H  = 


m=  1 
M 

©O^m  ®  Ip)BmS[t 

171=  1 


M 


©  (B™  +  ^Jp)  ©  (  (^m  0  Ip)BmSL 


771=1 


M 


771=1 

M 


771=1 


(J)  ©  ^p)  J  (J)  ( SLBmS £  +  ©j/ 


771=1 


Since  both  Bm  +  a^Ip  and  SLBmS l  +  cr^IL  are  positive  definite  matrices,  we  may  define: 


U  =  UT  = 


M 


771=1 


©(-Bm  +  (7^1  p)  1/2  0 

M 

0  Q}{SLBrnSr[  +  apL)~1/ 2 


771=1 


(B.3) 
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so  that 


U1  HU  = 


Imp 

m 

©Fri 

L  m= 1 


where 


M 


©  F” 
m=  1 

In 


Fm  =  (Bm  +  PnIp)-1'2(a^<DIp)BmSl(SLBmS[  +  PjL)-1'2. 


(B.4) 


Defining 


V  = 


M 


®F„ 


LMP  yjy  ^  m 

m= 1 

0  /a- 


(B.5) 


we  apply  Er  on  the  left  of  UT HU  and  V  on  the  right  to  obtain  the  block  diagonal  matrix 
D: 


D  =  (UV)TH(UV)  =  VT{UTHU)V 


Imp 


0 


0 

M 

®(4  -  F*F.n) 


(B.6) 


The  inverse  of  U  is  obtained  by  replacing  the  exponent  —1/2  with  1/2  in  Equation 
(B.3)  and  the  inverse  of  V  by  replacing  —  ©  with  ©  in  Equation  (B.5).  Therefore,  UV  is 
nonsingular  and  we  may  solve  for  H: 

H  =  [(UV)T]~l  D(UV)~l . 


We  see  that  H  is  nonsingular  if  and  only  if  D  is  nonsingular  and  in  that  case: 

H~l  =  ( UV)D~\UV)t . 

Thus,  determining  the  singular  or  nonsingular  status  of  H  reduces  to  determining  the  sin¬ 
gular  or  nonsingular  status  of  II  —  F^Fm  for  m  =  1,  2 M.  For  the  special  case  where 
the  P  x  P  diagonal  blocks  of  C~  are  identical  within  super-pixels,  this  allows  the  singular 
or  nonsingular  status  of  H  to  be  characterized  in  terms  of  an  and  av  as  indicated  in  the 
following  theorem: 


Theorem  B.0.5  For  the  special  case  Bmj  =  Bm^  (1  <j<  L,  1  <  m  <  M),  the  matrix  H 
of  Equation  (B.2)  is  positive  definite  (hence  non-singular)  when  either  an  >  0  or  av  >  0, 
and  is  singular  whenever  an  =  av  =  0. 
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Proof.  From  the  development  above,  H  is  positive  definite  if  and  only  if  A  —  FfnFm  is 
positive  definite  for  m  =  1,2,...,  M.  In  general, 


FZFm  =  (, 

<j= 1 


-*T 
(jJmUJ 

ai  m 


d= 1 


where  dm  j  =  sTBmjS  +  ai.  Under  the  special  case  assumption  this  reduces  to: 


T  jp  _  ®  1  “1“  ^n^P)  1® 


F  = 

1  m1  m 


STFm5iS  + 


(B-7) 


where  ym  =  cu^ci ;m  =  ||<AnA^||F-  Comparing  with  Equation  (A. 15),  we  see  that  F^Fm  =  Am. 

To  prove  that  A  —  FjnFrn  is  positive  definite,  let  u  be  any  unit  vector  of  length  L  and 
show 


ut(A  -  F^Fm) u  >  0. 

Since,  uT(A  —  FfLFm)u  =  1  —  ||Fmu|||,  this  reduces  to  showing  ||Fmu||2  <  1.  Using  the 
observation  that  the  2-norm  and  Frobenius  norm  are  identical  when  applied  to  vectors,  it 
follows  that: 

||AnU||2  <  ||Fm||i?  ||u||2  =  \\Fm\\F. 

By  Lemma  A. 0.1,  ||F^Fm||j?  <  1  and 

ll-^m-Cra||F  <  1  ==>■  WF^Wf  1 1  Fm  1 1  p  <  1  =>  \\Fj\f  <  1  =>  1 1  Fm  1 1  p  <  1. 


Therefore,  ||Tmu||  2  <  1-  ■ 

In  the  proof  of  Theorem  B.0.5,  it  is  shown  that  A  —  FffFm  =  IL  —  \m  is  positive  definite. 
It  follows  that  the  matrix  SRST  +  CV  appearing  in  the  Form  2  MAP  estimate  (see  Theorem 
A. 0.3)  not  only  is  nonsingular  under  the  special  case  assumption  and  the  assumption  that 
either  an  >  0  or  av  >  0,  it  is  in  fact  positive  definite  under  these  assumptions. 


We  now  derive  the  Form  3  MAP  estimate  for  the  noise  free  case  (an  =  av  =  0),  under 
the  special  case  assumption  that  the  P  x  P  blocks  of  Cf  are  identical  within  super-pixels. 
As  before,  we  take  an  =  0  and  examine  A(:i)  in  the  limit  as  crv  — »  0.  Thus  for  the  remainder 
of  this  section,  we  assume  an  =  0  and  Bmj  =  Am> i  for  1  <  m  <  M  and  1  <  j  <  L. 


We  successively  derive  expressions  for  A  1,VD  1VT ,  H  1  =  U(VD  lVT)UT  and  finally 


A®  = 


CzW1 ,  CzS1  H  ,  taking  the  limit  as  crr,  — »■  0  in  the  final  expression.  Using 


sL 

=  II  ®  sT, 

BmSl 

—  Il  &  -^m,  1^ 
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and  the  properties  of  the  Kronecker  tensor  product,  we  derive  from  Equation  (B.4)  the 
following  expression  for  Fm\ 


F  = 

±  m 


®  BmlS 


(B.8) 


where  dm  =  sTBm:iS  +  Either  of  Equations  (B.7)  or  (B.8)  yields: 

I L  d^mFrn  I L  Pm  i^rn^rn)  t 


where 


Pm 


dm  1 


Tm 


(as  (j^  — >  0). 


Defining  Em  =  IL-  we  have: 

M 


D~l  = 


0/p  o 

=i 

M 

o  ©iC 


m= 1 


m=l 


,  and 


VD~lVT  = 


Since  U  =  UT  is  given  by: 


M 


M 


0  (Ip  +  Fr*E-'FZ)  -©E^E"1 


771=1 


M 


771=1 

M 


©F,;'F,I  ©E-1 


771=1 


771=1 


U  = 


M 


©t y/2<r  o 


771=1 


M 


0  0  dm/2h 


771=1 


it  follows  that  H -1  =  U {V D~lVT)UT  is  given  by: 
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Using  Equation  (B.8)  and  properties  of  the  Kronecker  tensor  product,  this  reduces  to: 


H -1  = 


M 

©f 

m= 1 


B~\  + 


^mdf) 


)  (, 


UJm  ©  S  ) 


M 


-  © 

m= 1  7m“m 


Wra  ®  ST) 


“  ©  (^m  0  S)  ^m1 

m=i  7m  m 
M 


In  order  facilitate  the  matrix  algebra,  we  partition  A ®  into  blocks  of  size  NP  x  MP 


m= 1 

(3) 


and  of  size  iVP  x  N: 


A W  =  C~z 


W 

S 


H 


-l  _ 


4  (3)  4(3) 
Xll  )  Zl2 


Since 


cz- 


iu 

5 


M 


M 


0  ((Um  ©  Bm,l)  •>  ©(^  ®  Bmt is) 


.771—1 


771—1 
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Applying  the  properties  of  the  Kronecker  tensor  product,  these  reduce  to  the  following: 
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By  Lemma  A. 0.4: 


1 


0  Bm\S 

7m 


(  II  -  —  )  PmX  =  h-  —  (for  CT„  >  0). 

V  7 m  )  7m 


123 


Therefore,  for  av  >  0.  we  have: 
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Taking  the  limit  as  av  — >  0,  the  coefficient  matrix  A('3)  under  the  special  case  and  noise  free 
assumptions  is: 
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where  a,„  =  sTh>mls.  Breaking  Equation  (B.l)  into  super-pixel  components  yields: 
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where  is  the  component  of  Equation  (B.9).  Simplification  of  this  yields  the  rrfi1 
super-pixel  of  the  Form  3  MAP  estimate  under  the  special  case  and  noise  free  assumptions: 
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As  a  check  on  our  work,  we  see  that  Equations  (4.70),  (A. 23),  and  (B.10)  are  identical. 
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Appendix  C 


Temperature /Emissivity  Separation 
and  Atmospheric  Correction 

C.l  Introduction 

This  appendix  contains  a  literature  review  and  preliminary  experiments  done  to  compensate 
for  atmospheric  effects  in  hyperspectral  data.  This  preliminary  work  was  carried  out  by 
Precila  C.  Ip  at  Affiant  Techsystems  -  Mission  Research  in  Nashua,  New  Hampshire,  03062- 
1323.  It  is  not  intended  to  be  a  complete  study  of  atmospheric  effects,  but  rather  a  modest 
supplement  to  the  primary  effort  described  in  the  majority  of  the  report. 


C.2  Temperature/Emissivity  Separation  Algorithms 

C.2.1  Problem  Statement 

The  general  problem  of  Temperature/Emissivity  Separation  (TES)  is  a  challenging  one.  The 
challenge  arises  because  it  is  inherently  an  underdetermined  problem:  in  a  multispectral  or 
hyperspectral  radiance  measurement  with  N  channels,  there  are  N+l  unknowns.  There  is 
one  unknown  emissivity  per  channel  and  one  surface  temperature  unknown.  For  specific 
cases,  such  as  hyperspectral  data  collected  over  oceans  where  the  emissivity  of  water  is 
known,  the  problem  is  easily  solved.  However,  for  data  taken  over  land  surfaces  where  the 
emissivity  of  the  background  materials  are  unknown  or  highly  variable,  the  problem  is  solv¬ 
able  only  when  additional  assumptions  or  approximations  are  made  concerning  background 
materials  and  their  emissivity.  Many  iterative  solutions  have  been  proposed  that  constrain 
the  extra  degree  of  freedom. 
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C.2.2  Literature  Review 


In  1993,  Kealy  and  Hook  [29]  evaluated  three  methods  for  recovery  of  temperature  and 
emissivity  from  multispectral  thermal  infrared  (TIR)  radiance  spectra  acquired  over  land 
surfaces.  The  instruments  for  data  collection  for  this  study  were  the  Thermal  Infrared  Mul¬ 
tispectral  Scanner  (TIMS)  and  the  Advanced  Spaceborne  Thermal  Emission  Reflectance 
Radiometer  (ASTER).  ASTER  is  a  scanning  instrument  on  the  NASA’s  Terra  Satellite 
for  recording  thermal  IR  data  to  obtain  surface  temperatures  and  emissivity  spectra.  For 
ASTER,  there  are  5  channels,  therefore  the  radiance  estimation  problem  contains  5  measure¬ 
ments,  but  6  unknowns  when  surface  temperature  is  included.  The  three  methods  evaluated 
by  Kealy  and  Hook  were:  (1)  reference  channel,  (2)  normalized  emissivity  (NEM),  and  (3) 
alpha-derived  emissivity  (ADE).  The  result  of  their  study  showed  that  NEM  and  ADE  are 
more  accurate  than  the  reference  channel  method.  Between  NEM  and  ADE,  ADE  is  supe¬ 
rior  for  terrain  backgrounds  with  widely  varying  material  emissivity  such  as  vegetation  and 
igneous  rocks.  The  ADE  method  is  also  equally  accurate  for  data  that  was  convolved  with 
the  filter  response  functions  of  the  TIMS  and  ASTER  instruments. 

In  1994,  a  working  group  (ASTER  Temperature/Emissivity  Separation  Algorithm  Fea¬ 
sibility  Study)  was  formed  to  evaluate  ten  existing  algorithms  for  obtaining  temperatures 
and  emissivity.  For  each  of  the  ten  algorithms,  [30]  listed  the  reasons  for  rejection  of  the 
particular  algorithm  by  the  ASTER  working  group.  The  reviewed  algorithms  were: 

1.  Alpha-derived  emissivity  (ADE)  method 

2.  Classification  method 

3.  Day-night  measurement 

4.  Emissivity  bounds  method 

5.  Graybody  emissivity  method 

6.  Mean-MMD  method  (MMD) 

7.  Model  emissivity  method 

8.  Normalized  emissivity  method  (NEM) 

9.  Ratio  algorithm  (RAT) 

10.  Split-window  algorithm 

The  members  of  this  group  consisted  of  experts  in  the  field  from  U.S.  and  Japan.  Since 
this  initial  meeting,  the  group  has  met  regularly  to  develop,  improve,  and  validate  algorithms 
for  high  accuracy  land  surface  temperature  and  emissivity  computation.  There  other  pri¬ 
mary  function  is  to  distribute  the  products  to  the  user  community.  The  websites  for  the 
ASTER  data  products  and  the  working  group  are: 

http : //asterweb . jpl . nasa . gov/products/data_products . htm 

http :  //www .  science .  aster .  ersdac .  or .  jp/en/science_inf  o/te_WG .  html 
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Available  data  products  relevant  to  the  work  presented  here  are  AST08  (surface  kinetic 
temperature),  AST05  (surface  emissivity),  and  AST09T  (atmospheric  corrected  surface  ra¬ 
diance).  These  three  products  have  been  validated  for  cloud- free  pixels. 


C.2.3  The  Hybrid  TES  Algorithm 

Based  on  the  reviews,  the  ASTER  working  group  developed  a  hybrid  algorithm  based  on 
two  algorithms:  NEM  and  MMD  (the  latter  is  based  on  ADE).  The  advantages  of  this 
hybrid  TES  algorithm  over  NEM  or  MMD  alone  are  greater  accuracy  and  precision  for 
temperature  and  emissivity  recovery.  The  goals  of  the  hybrid  TES  algorithm  are  to  obtain 
surface  temperatures  especially  over  vegetation,  water,  and  snow  as  well  as  to  obtain  surface 
emissivity  for  mineral  substrates. 

In  2002,  Dash  et  al.  [31]  did  a  comprehensive  review  of  methods  on  estimating  land 
surface  temperature  (LST)  and  emissivity  using  passive  sensor  data.  He  concluded  that  the 
hybrid  TES  algorithm  is  applicable  and  viable  if  information  on  the  atmosphere  and  the 
multiple  IR  channels  are  available.  The  main  disadvantage  for  using  hybrid  TES  algorithm 
is  that  accurate  determination  of  the  downwelling  atmospheric  irradiance  is  critical  to  the 
results. 

Most  recently,  in  2004,  Payan  et  al.  [32]  evaluated  the  hybrid  TES  algorithm  for  both 
hyperspectral  and  multispectral  data.  They  confirmed  the  applicability  of  the  algorithm  in 
both  long-wave  infrared  (LWIR)  and  mid-wave  infrared  (MWIR)  for  shaded  surfaces.  How¬ 
ever,  this  study  [32]  found  that  the  methods  cannot  retrieve  emissivity  for  high  contrast 
surfaces  such  as  metal.  It  also  showed  that  the  algorithms  cannot  be  used  to  estimate  emis¬ 
sivity  for  hyperspectral  data  in  the  MWIR  for  cases  where  the  sun  is  directly  illuminating 
the  surface. 

Based  on  the  studies  and  review  in  [30-32] ,  we  concluded  that  the  hybrid  TES  algorithm 
is  the  most  robust  for  emissivity  and  temperature  calculation  in  hyperspectral  data.  In 
future  work,  we  would  like  to  implement  and  exhaustively  test  this  hybrid  algorithm  for 
various  temperature/emissivity  separation  tasks. 

The  steps  for  implementing  the  TES  algorithm,  which  is  a  combination  of  the  NEM, 
RAT,  and  MMD  modules  are  [30]: 

1.  Estimate  the  surface  temperature  and  subtract  the  reflected  sky  irradiance  iteratively 
(NEM).  This  method  first  removes  the  atmospheric  radiance  and  then  estimates  the 
temperature  and  emissivity  using  an  assumed  maximum  emissivity. 

2.  Obtain  relative  emissivities  and  spectral  shape  by  performing  a  ratio  of  the  emissivities 
from  NEM  to  the  average  emissivity  (RAT). 

3.  Estimate  the  TES  emissivities  and  temperature  using  regression  analysis  (MMD).  This 
method  obtains  the  minimum  emissivity  empirically  and  then  calculates  the  absolute 
emissivity. 
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Detailed  flow  charts  and  description  of  the  hybrid  TES  algorithm  is  available  in  [30]. 
This  algorithm  may  become  available  as  a  consumer  off-the-shelf  product  in  the  near  future. 


C.2.4  Preliminary  Temperature/Emissivity  Separation  Results 

As  a  preliminary  and  exploratory  experiment,  we  decided  to  verify  the  performance  of 
the  Reference  Channel  Method  and  the  Normalized  Emissivity  Method  used  by  Kealy  and 
Hook  [29]  on  a  representative  hyperspectral  data  cube.  To  perform  this  task,  we  used  ENVI 
(the  Environment  for  Visualizing  Images)  Version  3.6  to  generate  surface  temperatures  from 
radiance  data  and  apply  the  algorithms  described.  It  should  be  noted  that  ENVI  3.6  does 
not  generate  temperature  using  the  Alpha  Residual  algorithm  (the  third  method  discussed 
in  [29]),  although  it  is  likely  that  it  will  be  incorporated  into  a  future  release  of  ENVI. 

The  hyperspectral  data  set  utilized  was  taken  by  the  NASA  MODIS/ASTER  Airborne 
Simulator  (MASTER)  scanner  over  the  Cuprite  Mining  District  in  Nevada.  This  scanner 
covers  50  bands  in  the  spectral  range  0.46-12.845  microns.  For  retrieving  surface  temper¬ 
ature,  only  the  9  thermal  bands  from  8.2-12.845  microns  are  used.  The  results  of  these 
algorithms  are  shown  in  Figure  C.1-C.2.  Surface  temperature  (K)  statistics  derived  using 
these  algorithms  are  shown  in  Table  C.l. 


Figure  C.l:  The  temperature  histogram  for  a  716  by  2028  image  using  the  Reference  Channel 
algorithm  assuming  an  emissivity  of  0.96. 


Table  C.l:  Derived  Surface  Temperature  (K)  Statistics  versus  method  Algorithm 


Algorithm 

Min 

Max 

Mean 

Stdev 

Reference  Channel 

287.323212 

323.351776 

309.294547 

3.189591 

Normalized  Emissivity 

292.194824 

331.534668 

319.103284 

3.652123 
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Figure  C.2:  The  temperature  histogram  for  a  716x2028  image  using  the  Normalized  Emis- 
sivity  algorithm  assuming  an  emissivity  of  0.96. 

C.3  Atmospheric  Correction 

In  addition  to  investigating  Temperature/Emissivity  Separation  Algorithms,  we  explored  at¬ 
mospheric  correction  algorithms  that  might  be  integrated  into  the  MAP  algorithm  presented 
in  the  main  body  of  the  report.  Two  candidate  algorithms  were  identified:  FLAASH  (Fast 
Line-of-sight  Atmospheric  Analysis  of  Spectral  Hypercubes),  and  Parallel  MODTRAN.  The 
details  of  the  algorithms  are  not  presented  here,  but  are  available  in  [33]  and  [34], 

The  FLAASH  algorithm  may  be  used  to  retrieve  reflectance  data  from  radiance  data.  To 
demonstrate  this,  an  representative  example  of  AVIRIS  data  collected  over  Fort  AP  Hill  in 
Virginia  on  September  26,  2001,  was  processed  by  FLAASH  .  Figure  C.3  shows  the  radiance 
spectrum  from  400  to  2500  nm  for  a  vegetation  pixel,  and  Figure  C.4  shows  the  reflectance 
spectrum  after  FLAASH  was  applied  for  the  same  vegetation  pixel. 

For  future  work,  there  are  parallel  versions  of  FLAASH  and  MODTRAN  3  that  might 
be  incorporated  into  a  parallel  MAP  algorithm.  As  a  preliminary  step,  we  implemented  the 
NASA/JPL  parallel  version  of  MODTRAN  3  on  a  Beowulf  Cluster.  Figure  C.5  shows  the 
processing  time  versus  the  number  of  processors. 
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Figure  C.4:  AVIRIS  Reflectance  Spectrum  at  Fort  AP  Hill  obtained  from  FLAASH. 
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Figure  C.5:  Time  as  a  function  of  number  of  processors 


C.4  Summary 

In  this  appendix,  we  reviewed  and  identified  several  algorithms  for  temperature/emissivity 
separation.  The  hybrid  algorithm  described  in  Section  C.2.3  emerged  as  the  current  state-of- 
the-art  algorithm  for  this  purpose.  For  atmospheric  correction,  we  identified  parallel  versions 
of  the  FLAASH  and  MODTRAN  algorithms.  Future  work  would  include  incorporating 
parallel  implementations  of  the  hybrid  TES  algorithm  and  FLAASH,  as  parts  of  a  complete 
MAP  enhancement  of  hyperspectral  imagery. 
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